SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ pctrrfs()

subroutine pctrrfs ( character  uplo,
character  trans,
character  diag,
integer  n,
integer  nrhs,
complex, dimension( * )  a,
integer  ia,
integer  ja,
integer, dimension( * )  desca,
complex, dimension( * )  b,
integer  ib,
integer  jb,
integer, dimension( * )  descb,
complex, dimension( * )  x,
integer  ix,
integer  jx,
integer, dimension( * )  descx,
real, dimension( * )  ferr,
real, dimension( * )  berr,
complex, dimension( * )  work,
integer  lwork,
real, dimension( * )  rwork,
integer  lrwork,
integer  info 
)

Definition at line 1 of file pctrrfs.f.

4*
5* -- ScaLAPACK routine (version 1.7) --
6* University of Tennessee, Knoxville, Oak Ridge National Laboratory,
7* and University of California, Berkeley.
8* May 1, 1997
9*
10* .. Scalar Arguments ..
11 CHARACTER DIAG, TRANS, UPLO
12 INTEGER INFO, IA, IB, IX, JA, JB, JX, LRWORK, LWORK,
13 $ N, NRHS
14* ..
15* .. Array Arguments ..
16 INTEGER DESCA( * ), DESCB( * ), DESCX( * )
17 REAL BERR( * ), FERR( * ), RWORK( * )
18 COMPLEX A( * ), B( * ), WORK( * ), X( * )
19* ..
20*
21* Purpose
22* =======
23*
24* PCTRRFS provides error bounds and backward error estimates for the
25* solution to a system of linear equations with a triangular
26* coefficient matrix.
27*
28* The solution matrix X must be computed by PCTRTRS or some other
29* means before entering this routine. PCTRRFS does not do iterative
30* refinement because doing so cannot improve the backward error.
31*
32* Notes
33* =====
34*
35* Each global data object is described by an associated description
36* vector. This vector stores the information required to establish
37* the mapping between an object element and its corresponding process
38* and memory location.
39*
40* Let A be a generic term for any 2D block cyclicly distributed array.
41* Such a global array has an associated description vector DESCA.
42* In the following comments, the character _ should be read as
43* "of the global array".
44*
45* NOTATION STORED IN EXPLANATION
46* --------------- -------------- --------------------------------------
47* DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
48* DTYPE_A = 1.
49* CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
50* the BLACS process grid A is distribu-
51* ted over. The context itself is glo-
52* bal, but the handle (the integer
53* value) may vary.
54* M_A (global) DESCA( M_ ) The number of rows in the global
55* array A.
56* N_A (global) DESCA( N_ ) The number of columns in the global
57* array A.
58* MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
59* the rows of the array.
60* NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
61* the columns of the array.
62* RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
63* row of the array A is distributed.
64* CSRC_A (global) DESCA( CSRC_ ) The process column over which the
65* first column of the array A is
66* distributed.
67* LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
68* array. LLD_A >= MAX(1,LOCr(M_A)).
69*
70* Let K be the number of rows or columns of a distributed matrix,
71* and assume that its process grid has dimension p x q.
72* LOCr( K ) denotes the number of elements of K that a process
73* would receive if K were distributed over the p processes of its
74* process column.
75* Similarly, LOCc( K ) denotes the number of elements of K that a
76* process would receive if K were distributed over the q processes of
77* its process row.
78* The values of LOCr() and LOCc() may be determined via a call to the
79* ScaLAPACK tool function, NUMROC:
80* LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
81* LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
82* An upper bound for these quantities may be computed by:
83* LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
84* LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
85*
86* In the following comments, sub( A ), sub( X ) and sub( B ) denote
87* respectively A(IA:IA+N-1,JA:JA+N-1), X(IX:IX+N-1,JX:JX+NRHS-1) and
88* B(IB:IB+N-1,JB:JB+NRHS-1).
89*
90* Arguments
91* =========
92*
93* UPLO (global input) CHARACTER*1
94* = 'U': sub( A ) is upper triangular;
95* = 'L': sub( A ) is lower triangular.
96*
97* TRANS (global input) CHARACTER*1
98* Specifies the form of the system of equations.
99* = 'N': sub( A ) * sub( X ) = sub( B ) (No transpose)
100* = 'T': sub( A )**T * sub( X ) = sub( B ) (Transpose)
101* = 'C': sub( A )**H * sub( X ) = sub( B )
102* (Conjugate transpose)
103*
104* DIAG (global input) CHARACTER*1
105* = 'N': sub( A ) is non-unit triangular;
106* = 'U': sub( A ) is unit triangular.
107*
108* N (global input) INTEGER
109* The order of the matrix sub( A ). N >= 0.
110*
111* NRHS (global input) INTEGER
112* The number of right hand sides, i.e., the number of columns
113* of the matrices sub( B ) and sub( X ). NRHS >= 0.
114*
115* A (local input) COMPLEX pointer into the local memory
116* to an array of local dimension (LLD_A,LOCc(JA+N-1) ). This
117* array contains the local pieces of the original triangular
118* distributed matrix sub( A ).
119* If UPLO = 'U', the leading N-by-N upper triangular part of
120* sub( A ) contains the upper triangular part of the matrix,
121* and its strictly lower triangular part is not referenced.
122* If UPLO = 'L', the leading N-by-N lower triangular part of
123* sub( A ) contains the lower triangular part of the distribu-
124* ted matrix, and its strictly upper triangular part is not
125* referenced.
126* If DIAG = 'U', the diagonal elements of sub( A ) are also
127* not referenced and are assumed to be 1.
128*
129* IA (global input) INTEGER
130* The row index in the global array A indicating the first
131* row of sub( A ).
132*
133* JA (global input) INTEGER
134* The column index in the global array A indicating the
135* first column of sub( A ).
136*
137* DESCA (global and local input) INTEGER array of dimension DLEN_.
138* The array descriptor for the distributed matrix A.
139*
140* B (local input) COMPLEX pointer into the local memory
141* to an array of local dimension (LLD_B, LOCc(JB+NRHS-1) ).
142* On entry, this array contains the the local pieces of the
143* right hand sides sub( B ).
144*
145* IB (global input) INTEGER
146* The row index in the global array B indicating the first
147* row of sub( B ).
148*
149* JB (global input) INTEGER
150* The column index in the global array B indicating the
151* first column of sub( B ).
152*
153* DESCB (global and local input) INTEGER array of dimension DLEN_.
154* The array descriptor for the distributed matrix B.
155*
156* X (local input) COMPLEX pointer into the local memory
157* to an array of local dimension (LLD_X, LOCc(JX+NRHS-1) ).
158* On entry, this array contains the the local pieces of the
159* solution vectors sub( X ).
160*
161* IX (global input) INTEGER
162* The row index in the global array X indicating the first
163* row of sub( X ).
164*
165* JX (global input) INTEGER
166* The column index in the global array X indicating the
167* first column of sub( X ).
168*
169* DESCX (global and local input) INTEGER array of dimension DLEN_.
170* The array descriptor for the distributed matrix X.
171*
172* FERR (local output) REAL array of local dimension
173* LOCc(JB+NRHS-1). The estimated forward error bounds for
174* each solution vector of sub( X ). If XTRUE is the true
175* solution, FERR bounds the magnitude of the largest entry
176* in (sub( X ) - XTRUE) divided by the magnitude of the
177* largest entry in sub( X ). The estimate is as reliable as
178* the estimate for RCOND, and is almost always a slight
179* overestimate of the true error.
180* This array is tied to the distributed matrix X.
181*
182* BERR (local output) REAL array of local dimension
183* LOCc(JB+NRHS-1). The componentwise relative backward
184* error of each solution vector (i.e., the smallest re-
185* lative change in any entry of sub( A ) or sub( B )
186* that makes sub( X ) an exact solution).
187* This array is tied to the distributed matrix X.
188*
189* WORK (local workspace/local output) COMPLEX array,
190* dimension (LWORK)
191* On exit, WORK(1) returns the minimal and optimal LWORK.
192*
193* LWORK (local or global input) INTEGER
194* The dimension of the array WORK.
195* LWORK is local input and must be at least
196* LWORK >= 2*LOCr( N + MOD( IA-1, MB_A ) ).
197*
198* If LWORK = -1, then LWORK is global input and a workspace
199* query is assumed; the routine only calculates the minimum
200* and optimal size for all work arrays. Each of these
201* values is returned in the first entry of the corresponding
202* work array, and no error message is issued by PXERBLA.
203*
204* RWORK (local workspace/local output) REAL array,
205* dimension (LRWORK)
206* On exit, RWORK(1) returns the minimal and optimal LRWORK.
207*
208* LRWORK (local or global input) INTEGER
209* The dimension of the array RWORK.
210* LRWORK is local input and must be at least
211* LRWORK >= LOCr( N + MOD( IB-1, MB_B ) ).
212*
213* If LRWORK = -1, then LRWORK is global input and a workspace
214* query is assumed; the routine only calculates the minimum
215* and optimal size for all work arrays. Each of these
216* values is returned in the first entry of the corresponding
217* work array, and no error message is issued by PXERBLA.
218*
219*
220* INFO (global output) INTEGER
221* = 0: successful exit
222* < 0: If the i-th argument is an array and the j-entry had
223* an illegal value, then INFO = -(i*100+j), if the i-th
224* argument is a scalar and had an illegal value, then
225* INFO = -i.
226*
227* Notes
228* =====
229*
230* This routine temporarily returns when N <= 1.
231*
232* The distributed submatrices sub( X ) and sub( B ) should be
233* distributed the same way on the same processes. These conditions
234* ensure that sub( X ) and sub( B ) are "perfectly" aligned.
235*
236* Moreover, this routine requires the distributed submatrices sub( A ),
237* sub( X ), and sub( B ) to be aligned on a block boundary,
238* i.e., if f(x,y) = MOD( x-1, y ):
239* f( IA, DESCA( MB_ ) ) = f( JA, DESCA( NB_ ) ) = 0,
240* f( IB, DESCB( MB_ ) ) = f( JB, DESCB( NB_ ) ) = 0, and
241* f( IX, DESCX( MB_ ) ) = f( JX, DESCX( NB_ ) ) = 0.
242*
243* =====================================================================
244*
245* .. Parameters ..
246 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
247 $ LLD_, MB_, M_, NB_, N_, RSRC_
248 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
249 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
250 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
251 REAL ZERO, RONE
252 parameter( zero = 0.0e+0, rone = 1.0e+0 )
253 COMPLEX ONE
254 parameter( one = ( 1.0e+0, 0.0e+0 ) )
255* ..
256* .. Local Scalars ..
257 LOGICAL LQUERY, NOTRAN, NOUNIT, UPPER
258 CHARACTER TRANSN, TRANST
259 INTEGER IAROW, IXBCOL, IXBROW, IXCOL, IXROW, ICOFFA,
260 $ ICOFFB, ICOFFX, ICTXT, ICURCOL, IDUM, II, IIXB,
261 $ IIW, IOFFXB, IPB, IPR, IPV, IROFFA, IROFFB,
262 $ IROFFX, IW, J, JBRHS, JJ, JJFBE, JJXB, JN, JW,
263 $ K, KASE, LDXB, LRWMIN, LWMIN, MYCOL, MYRHS,
264 $ MYROW, NP, NP0, NPCOL, NPMOD, NPROW, NZ
265 REAL EPS, EST, LSTRES, S, SAFE1, SAFE2, SAFMIN
266 COMPLEX ZDUM
267* ..
268* .. Local Arrays ..
269 INTEGER DESCW( DLEN_ ), IDUM1( 5 ), IDUM2( 5 )
270* ..
271* .. External Functions ..
272 LOGICAL LSAME
273 INTEGER ICEIL, INDXG2P, NUMROC
274 REAL PSLAMCH
275 EXTERNAL iceil, indxg2p, lsame, numroc, pslamch
276* ..
277* .. External Subroutines ..
278 EXTERNAL blacs_gridinfo, cgebr2d, cgebs2d, chk1mat,
279 $ descset, infog2l, pcatrmv, pcaxpy, pchk1mat,
280 $ pchk2mat, pccopy, pclacon, pctrmv,
281 $ pctrsv, pxerbla, sgamx2d
282* ..
283* .. Intrinsic Functions ..
284 INTRINSIC abs, aimag, cmplx, ichar, max, min, mod, real
285* ..
286* .. Statement Functions ..
287 REAL CABS1
288* ..
289* .. Statement Function definitions ..
290 cabs1( zdum ) = abs( real( zdum ) ) + abs( aimag( zdum ) )
291* ..
292* .. Executable Statements ..
293*
294* Get grid parameters
295*
296 ictxt = desca( ctxt_ )
297 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
298*
299* Test the input parameters.
300*
301 info = 0
302 IF( nprow.EQ.-1 ) THEN
303 info = -( 900+ctxt_ )
304 ELSE
305 CALL chk1mat( n, 4, n, 4, ia, ja, desca, 9, info )
306 CALL chk1mat( n, 4, nrhs, 5, ib, jb, descb, 13, info )
307 CALL chk1mat( n, 4, nrhs, 5, ix, jx, descx, 17, info )
308 IF( info.EQ.0 ) THEN
309 upper = lsame( uplo, 'U' )
310 notran = lsame( trans, 'N' )
311 nounit = lsame( diag, 'N' )
312 iroffa = mod( ia-1, desca( mb_ ) )
313 icoffa = mod( ja-1, desca( nb_ ) )
314 iroffb = mod( ib-1, descb( mb_ ) )
315 icoffb = mod( jb-1, descb( nb_ ) )
316 iroffx = mod( ix-1, descx( mb_ ) )
317 icoffx = mod( jx-1, descx( nb_ ) )
318 iarow = indxg2p( ia, desca( mb_ ), myrow, desca( rsrc_ ),
319 $ nprow )
320 CALL infog2l( ib, jb, descb, nprow, npcol, myrow, mycol,
321 $ iixb, jjxb, ixbrow, ixbcol )
322 ixrow = indxg2p( ix, descx( mb_ ), myrow, descx( rsrc_ ),
323 $ nprow )
324 ixcol = indxg2p( jx, descx( nb_ ), mycol, descx( csrc_ ),
325 $ npcol )
326 npmod = numroc( n+iroffa, desca( mb_ ), myrow, iarow,
327 $ nprow )
328 lwmin = 2*npmod
329 work( 1 ) = real( lwmin )
330 lrwmin = npmod
331 rwork( 1 ) = real( lrwmin )
332 lquery = ( lwork.EQ.-1 .OR. lrwork.EQ.-1 )
333*
334 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
335 info = -1
336 ELSE IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND.
337 $ .NOT.lsame( trans, 'C' ) ) THEN
338 info = -2
339 ELSE IF( .NOT.nounit .AND. .NOT.lsame( diag, 'U' ) ) THEN
340 info = -3
341 ELSE IF( n.LT.0 ) THEN
342 info = -4
343 ELSE IF( nrhs.LT.0 ) THEN
344 info = -5
345 ELSE IF( iroffa.NE.0 ) THEN
346 info = -7
347 ELSE IF( icoffa.NE.0 ) THEN
348 info = -8
349 ELSE IF( desca( mb_ ).NE.desca( nb_ ) ) THEN
350 info = -( 900+nb_ )
351 ELSE IF( iroffa.NE.iroffb .OR. iarow.NE.ixbrow ) THEN
352 info = -11
353 ELSE IF( desca( mb_ ).NE.descb( mb_ ) ) THEN
354 info = -( 1300+mb_ )
355 ELSE IF( ictxt.NE.descb( ctxt_ ) ) THEN
356 info = -( 1300+ctxt_ )
357 ELSE IF( iroffx.NE.0 .OR. ixbrow.NE.ixrow ) THEN
358 info = -15
359 ELSE IF( icoffb.NE.icoffx .OR. ixbcol.NE.ixcol ) THEN
360 info = -16
361 ELSE IF( descb( mb_ ).NE.descx( mb_ ) ) THEN
362 info = -( 1700+mb_ )
363 ELSE IF( descb( nb_ ).NE.descx( nb_ ) ) THEN
364 info = -( 1700+nb_ )
365 ELSE IF( ictxt.NE.descx( ctxt_ ) ) THEN
366 info = -( 1700+ctxt_ )
367 ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
368 info = -21
369 ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
370 info = -23
371 END IF
372 END IF
373*
374 IF( upper ) THEN
375 idum1( 1 ) = ichar( 'U' )
376 ELSE
377 idum1( 1 ) = ichar( 'L' )
378 END IF
379 idum2( 1 ) = 1
380 IF( notran ) THEN
381 idum1( 2 ) = ichar( 'N' )
382 ELSE IF( lsame( trans, 'T' ) ) THEN
383 idum1( 2 ) = ichar( 'T' )
384 ELSE
385 idum1( 2 ) = ichar( 'C' )
386 END IF
387 idum2( 2 ) = 2
388 IF( nounit ) THEN
389 idum1( 3 ) = ichar( 'N' )
390 ELSE
391 idum1( 3 ) = ichar( 'U' )
392 END IF
393 idum2( 3 ) = 3
394 IF( lwork.EQ.-1 ) THEN
395 idum1( 4 ) = -1
396 ELSE
397 idum1( 4 ) = 1
398 END IF
399 idum2( 4 ) = 21
400 IF( lrwork.EQ.-1 ) THEN
401 idum1( 5 ) = -1
402 ELSE
403 idum1( 5 ) = 1
404 END IF
405 idum2( 5 ) = 23
406 CALL pchk1mat( n, 4, n, 4, ia, ja, desca, 9, 0, idum1, idum2,
407 $ info )
408 CALL pchk2mat( n, 4, nrhs, 5, ib, jb, descb, 13, n, 4, nrhs, 5,
409 $ ix, jx, descx, 17, 5, idum1, idum2, info )
410 END IF
411 IF( info.NE.0 ) THEN
412 CALL pxerbla( ictxt, 'PCTRRFS', -info )
413 RETURN
414 ELSE IF( lquery ) THEN
415 RETURN
416 END IF
417*
418 jjfbe = jjxb
419 myrhs = numroc( jb+nrhs-1, descb( nb_ ), mycol, descb( csrc_ ),
420 $ npcol )
421*
422* Quick return if possible
423*
424 IF( n.LE.1 .OR. nrhs.EQ.0 ) THEN
425 DO 10 jj = jjfbe, myrhs
426 ferr( jj ) = zero
427 berr( jj ) = zero
428 10 CONTINUE
429 RETURN
430 END IF
431*
432 IF( notran ) THEN
433 transn = 'N'
434 transt = 'C'
435 ELSE
436 transn = 'C'
437 transt = 'N'
438 END IF
439*
440 np0 = numroc( n+iroffb, descb( mb_ ), myrow, ixbrow, nprow )
441 CALL descset( descw, n+iroffb, 1, desca( mb_ ), 1, ixbrow, ixbcol,
442 $ ictxt, max( 1, np0 ) )
443 ipb = 1
444 ipr = 1
445 ipv = ipr + np0
446 IF( myrow.EQ.ixbrow ) THEN
447 iiw = 1 + iroffb
448 np = np0 - iroffb
449 ELSE
450 iiw = 1
451 np = np0
452 END IF
453 iw = 1 + iroffb
454 jw = 1
455 ldxb = descb( lld_ )
456 ioffxb = ( jjxb-1 )*ldxb
457*
458* NZ = maximum number of nonzero entries in each row of A, plus 1
459*
460 nz = n + 1
461 eps = pslamch( ictxt, 'Epsilon' )
462 safmin = pslamch( ictxt, 'Safe minimum' )
463 safe1 = nz*safmin
464 safe2 = safe1 / eps
465 jn = min( iceil( jb, descb( nb_ ) )*descb( nb_ ), jb+nrhs-1 )
466*
467* Handle first block separately
468*
469 jbrhs = jn - jb + 1
470 DO 90 k = 0, jbrhs - 1
471*
472* Compute residual R = B - op(A) * X,
473* where op(A) = A or A', depending on TRANS.
474*
475 CALL pccopy( n, x, ix, jx+k, descx, 1, work( ipr ), iw, jw,
476 $ descw, 1 )
477 CALL pctrmv( uplo, trans, diag, n, a, ia, ja, desca,
478 $ work( ipr ), iw, jw, descw, 1 )
479 CALL pcaxpy( n, -one, b, ib, jb+k, descb, 1, work( ipr ), iw,
480 $ jw, descw, 1 )
481*
482* Compute componentwise relative backward error from formula
483*
484* max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
485*
486* where abs(Z) is the componentwise absolute value of the matrix
487* or vector Z. If the i-th component of the denominator is less
488* than SAFE2, then SAFE1 is added to the i-th components of the
489* numerator and denominator before dividing.
490*
491 IF( mycol.EQ.ixbcol ) THEN
492 IF( np.GT.0 ) THEN
493 DO 20 ii = iixb, iixb + np - 1
494 rwork( iiw+ii-iixb ) = cabs1( b( ii+ioffxb ) )
495 20 CONTINUE
496 END IF
497 END IF
498*
499 CALL pcatrmv( uplo, trans, diag, n, rone, a, ia, ja, desca, x,
500 $ ix, jx+k, descx, 1, rone, rwork( ipb ), iw, jw,
501 $ descw, 1 )
502*
503 s = zero
504 IF( mycol.EQ.ixbcol ) THEN
505 IF( np.GT.0 ) THEN
506 DO 30 ii = iiw - 1, iiw + np - 2
507 IF( rwork( ipb+ii ).GT.safe2 ) THEN
508 s = max( s, cabs1( work( ipr+ii ) ) /
509 $ rwork( ipb+ii ) )
510 ELSE
511 s = max( s, ( cabs1( work( ipr+ii ) )+safe1 ) /
512 $ ( rwork( ipb+ii )+safe1 ) )
513 END IF
514 30 CONTINUE
515 END IF
516 END IF
517*
518 CALL sgamx2d( ictxt, 'All', ' ', 1, 1, s, 1, idum, idum, 1,
519 $ -1, mycol )
520 IF( mycol.EQ.ixbcol )
521 $ berr( jjfbe ) = s
522*
523* Bound error from formula
524*
525* norm(X - XTRUE) / norm(X) .le. FERR =
526* norm( abs(inv(op(A)))*
527* ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
528*
529* where
530* norm(Z) is the magnitude of the largest component of Z
531* inv(op(A)) is the inverse of op(A)
532* abs(Z) is the componentwise absolute value of the matrix or
533* vector Z
534* NZ is the maximum number of nonzeros in any row of A, plus 1
535* EPS is machine epsilon
536*
537* The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
538* is incremented by SAFE1 if the i-th component of
539* abs(op(A))*abs(X) + abs(B) is less than SAFE2.
540*
541* Use PCLACON to estimate the infinity-norm of the matrix
542* inv(op(A)) * diag(W),
543* where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
544*
545 IF( mycol.EQ.ixbcol ) THEN
546 IF( np.GT.0 ) THEN
547 DO 40 ii = iiw - 1, iiw + np - 2
548 IF( rwork( ipb+ii ).GT.safe2 ) THEN
549 rwork( ipb+ii ) = cabs1( work( ipr+ii ) ) +
550 $ nz*eps*rwork( ipb+ii )
551 ELSE
552 rwork( ipb+ii ) = cabs1( work( ipr+ii ) ) +
553 $ nz*eps*rwork( ipb+ii ) + safe1
554 END IF
555 40 CONTINUE
556 END IF
557 END IF
558*
559 kase = 0
560 50 CONTINUE
561 IF( mycol.EQ.ixbcol ) THEN
562 CALL cgebs2d( ictxt, 'Rowwise', ' ', np, 1, work( ipr ),
563 $ descw( lld_ ) )
564 ELSE
565 CALL cgebr2d( ictxt, 'Rowwise', ' ', np, 1, work( ipr ),
566 $ descw( lld_ ), myrow, ixbcol )
567 END IF
568 descw( csrc_ ) = mycol
569 CALL pclacon( n, work( ipv ), iw, jw, descw, work( ipr ),
570 $ iw, jw, descw, est, kase )
571 descw( csrc_ ) = ixbcol
572*
573 IF( kase.NE.0 ) THEN
574 IF( kase.EQ.1 ) THEN
575*
576* Multiply by diag(W)*inv(op(A)').
577*
578 CALL pctrsv( uplo, transt, diag, n, a, ia, ja, desca,
579 $ work( ipr ), iw, jw, descw, 1 )
580 IF( mycol.EQ.ixbcol ) THEN
581 IF( np.GT.0 ) THEN
582 DO 60 ii = iiw - 1, iiw + np - 2
583 work( ipr+ii ) = rwork( ipb+ii )*work( ipr+ii )
584 60 CONTINUE
585 END IF
586 END IF
587 ELSE
588*
589* Multiply by inv(op(A))*diag(W).
590*
591 IF( mycol.EQ.ixbcol ) THEN
592 IF( np.GT.0 ) THEN
593 DO 70 ii = iiw - 1, iiw + np - 2
594 work( ipr+ii ) = rwork( ipb+ii )*work( ipr+ii )
595 70 CONTINUE
596 END IF
597 END IF
598 CALL pctrsv( uplo, transn, diag, n, a, ia, ja, desca,
599 $ work( ipr ), iw, jw, descw, 1 )
600 END IF
601 GO TO 50
602 END IF
603*
604* Normalize error.
605*
606 lstres = zero
607 IF( mycol.EQ.ixbcol ) THEN
608 IF( np.GT.0 ) THEN
609 DO 80 ii = iixb, iixb + np - 1
610 lstres = max( lstres, cabs1( x( ioffxb+ii ) ) )
611 80 CONTINUE
612 END IF
613 CALL sgamx2d( ictxt, 'Column', ' ', 1, 1, lstres, 1, idum,
614 $ idum, 1, -1, mycol )
615 IF( lstres.NE.zero )
616 $ ferr( jjfbe ) = est / lstres
617*
618 jjxb = jjxb + 1
619 jjfbe = jjfbe + 1
620 ioffxb = ioffxb + ldxb
621*
622 END IF
623*
624 90 CONTINUE
625*
626 icurcol = mod( ixbcol+1, npcol )
627*
628* Do for each right hand side
629*
630 DO 180 j = jn + 1, jb + nrhs - 1, descb( nb_ )
631 jbrhs = min( jb+nrhs-j, descb( nb_ ) )
632 descw( csrc_ ) = icurcol
633*
634 DO 170 k = 0, jbrhs - 1
635*
636* Compute residual R = B - op(A) * X,
637* where op(A) = A or A', depending on TRANS.
638*
639 CALL pccopy( n, x, ix, j+k, descx, 1, work( ipr ), iw, jw,
640 $ descw, 1 )
641 CALL pctrmv( uplo, trans, diag, n, a, ia, ja, desca,
642 $ work( ipr ), iw, jw, descw, 1 )
643 CALL pcaxpy( n, -one, b, ib, j+k, descb, 1, work( ipr ),
644 $ iw, jw, descw, 1 )
645*
646* Compute componentwise relative backward error from formula
647*
648* max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
649*
650* where abs(Z) is the componentwise absolute value of the
651* matrix or vector Z. If the i-th component of the
652* denominator is less than SAFE2, then SAFE1 is added to the
653* i-th components of the numerator and denominator before
654* dividing.
655*
656 IF( mycol.EQ.ixbcol ) THEN
657 IF( np.GT.0 ) THEN
658 DO 100 ii = iixb, iixb + np - 1
659 rwork( iiw+ii-iixb ) = cabs1( b( ii+ioffxb ) )
660 100 CONTINUE
661 END IF
662 END IF
663*
664 CALL pcatrmv( uplo, trans, diag, n, rone, a, ia, ja, desca,
665 $ x, ix, j+k, descx, 1, rone, rwork( ipb ), iw,
666 $ jw, descw, 1 )
667*
668 s = zero
669 IF( mycol.EQ.ixbcol ) THEN
670 IF( np.GT.0 ) THEN
671 DO 110 ii = iiw - 1, iiw + np - 2
672 IF( rwork( ipb+ii ).GT.safe2 ) THEN
673 s = max( s, cabs1( work( ipr+ii ) ) /
674 $ rwork( ipb+ii ) )
675 ELSE
676 s = max( s, ( cabs1( work( ipr+ii ) )+safe1 ) /
677 $ ( rwork( ipb+ii )+safe1 ) )
678 END IF
679 110 CONTINUE
680 END IF
681 END IF
682*
683 CALL sgamx2d( ictxt, 'All', ' ', 1, 1, s, 1, idum, idum, 1,
684 $ -1, mycol )
685 IF( mycol.EQ.ixbcol )
686 $ berr( jjfbe ) = s
687*
688* Bound error from formula
689*
690* norm(X - XTRUE) / norm(X) .le. FERR =
691* norm( abs(inv(op(A)))*
692* ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))/norm(X)
693*
694* where
695* norm(Z) is the magnitude of the largest component of Z
696* inv(op(A)) is the inverse of op(A)
697* abs(Z) is the componentwise absolute value of the matrix
698* or vector Z
699* NZ is the maximum number of nonzeros in any row of A,
700* plus 1
701* EPS is machine epsilon
702*
703* The i-th component of
704* abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
705* is incremented by SAFE1 if the i-th component of
706* abs(op(A))*abs(X) + abs(B) is less than SAFE2.
707*
708* Use PCLACON to estimate the infinity-norm of the matrix
709* inv(op(A)) * diag(W),
710* where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
711*
712 IF( mycol.EQ.ixbcol ) THEN
713 IF( np.GT.0 ) THEN
714 DO 120 ii = iiw - 1, iiw + np - 2
715 IF( rwork( ipb+ii ).GT.safe2 ) THEN
716 rwork( ipb+ii ) = cabs1( work( ipr+ii ) ) +
717 $ nz*eps*rwork( ipb+ii )
718 ELSE
719 rwork( ipb+ii ) = cabs1( work( ipr+ii ) ) +
720 $ nz*eps*rwork( ipb+ii ) + safe1
721 END IF
722 120 CONTINUE
723 END IF
724 END IF
725*
726 kase = 0
727 130 CONTINUE
728 IF( mycol.EQ.ixbcol ) THEN
729 CALL cgebs2d( ictxt, 'Rowwise', ' ', np, 1, work( ipr ),
730 $ descw( lld_ ) )
731 ELSE
732 CALL cgebr2d( ictxt, 'Rowwise', ' ', np, 1, work( ipr ),
733 $ descw( lld_ ), myrow, ixbcol )
734 END IF
735 descw( csrc_ ) = mycol
736 CALL pclacon( n, work( ipv ), iw, jw, descw, work( ipr ),
737 $ iw, jw, descw, est, kase )
738 descw( csrc_ ) = ixbcol
739*
740 IF( kase.NE.0 ) THEN
741 IF( kase.EQ.1 ) THEN
742*
743* Multiply by diag(W)*inv(op(A)').
744*
745 CALL pctrsv( uplo, transt, diag, n, a, ia, ja, desca,
746 $ work( ipr ), iw, jw, descw, 1 )
747 IF( mycol.EQ.ixbcol ) THEN
748 IF( np.GT.0 ) THEN
749 DO 140 ii = iiw - 1, iiw + np - 2
750 work( ipr+ii ) = rwork( ipb+ii )*
751 $ work( ipr+ii )
752 140 CONTINUE
753 END IF
754 END IF
755 ELSE
756*
757* Multiply by inv(op(A))*diag(W).
758*
759 IF( mycol.EQ.ixbcol ) THEN
760 IF( np.GT.0 ) THEN
761 DO 150 ii = iiw - 1, iiw + np - 2
762 work( ipr+ii ) = rwork( ipb+ii )*
763 $ work( ipr+ii )
764 150 CONTINUE
765 END IF
766 END IF
767 CALL pctrsv( uplo, transn, diag, n, a, ia, ja, desca,
768 $ work( ipr ), iw, jw, descw, 1 )
769 END IF
770 GO TO 130
771 END IF
772*
773* Normalize error.
774*
775 lstres = zero
776 IF( mycol.EQ.ixbcol ) THEN
777 IF( np.GT.0 ) THEN
778 DO 160 ii = iixb, iixb + np - 1
779 lstres = max( lstres, cabs1( x( ioffxb+ii ) ) )
780 160 CONTINUE
781 END IF
782 CALL sgamx2d( ictxt, 'Column', ' ', 1, 1, lstres, 1,
783 $ idum, idum, 1, -1, mycol )
784 IF( lstres.NE.zero )
785 $ ferr( jjfbe ) = est / lstres
786*
787 jjxb = jjxb + 1
788 jjfbe = jjfbe + 1
789 ioffxb = ioffxb + ldxb
790*
791 END IF
792*
793 170 CONTINUE
794*
795 icurcol = mod( icurcol+1, npcol )
796*
797 180 CONTINUE
798*
799 work( 1 ) = cmplx( real( lwmin ) )
800 rwork( 1 ) = real( lrwmin )
801*
802 RETURN
803*
804* End of PCTRRFS
805*
cmplx
float cmplx[2]
Definition pblas.h:136
chk1mat
subroutine chk1mat(ma, mapos0, na, napos0, ia, ja, desca, descapos0, info)
Definition chk1mat.f:3
descset
subroutine descset(desc, m, n, mb, nb, irsrc, icsrc, ictxt, lld)
Definition descset.f:3
iceil
integer function iceil(inum, idenom)
Definition iceil.f:2
indxg2p
integer function indxg2p(indxglob, nb, iproc, isrcproc, nprocs)
Definition indxg2p.f:2
infog2l
subroutine infog2l(grindx, gcindx, desc, nprow, npcol, myrow, mycol, lrindx, lcindx, rsrc, csrc)
Definition infog2l.f:3
numroc
integer function numroc(n, nb, iproc, isrcproc, nprocs)
Definition numroc.f:2
pslamch
real function pslamch(ictxt, cmach)
Definition pcblastst.f:7455
max
#define max(A, B)
Definition pcgemr.c:180
min
#define min(A, B)
Definition pcgemr.c:181
pchk1mat
subroutine pchk1mat(ma, mapos0, na, napos0, ia, ja, desca, descapos0, nextra, ex, expos, info)
Definition pchkxmat.f:3
pchk2mat
subroutine pchk2mat(ma, mapos0, na, napos0, ia, ja, desca, descapos0, mb, mbpos0, nb, nbpos0, ib, jb, descb, descbpos0, nextra, ex, expos, info)
Definition pchkxmat.f:175
pclacon
subroutine pclacon(n, v, iv, jv, descv, x, ix, jx, descx, est, kase)
Definition pclacon.f:3
pxerbla
subroutine pxerbla(ictxt, srname, info)
Definition pxerbla.f:2
lsame
logical function lsame(ca, cb)
Definition tools.f:1724
Here is the call graph for this function: