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

◆ pdtrrfs()

subroutine pdtrrfs ( character  uplo,
character  trans,
character  diag,
integer  n,
integer  nrhs,
double precision, dimension( * )  a,
integer  ia,
integer  ja,
integer, dimension( * )  desca,
double precision, dimension( * )  b,
integer  ib,
integer  jb,
integer, dimension( * )  descb,
double precision, dimension( * )  x,
integer  ix,
integer  jx,
integer, dimension( * )  descx,
double precision, dimension( * )  ferr,
double precision, dimension( * )  berr,
double precision, dimension( * )  work,
integer  lwork,
integer, dimension( * )  iwork,
integer  liwork,
integer  info 
)

Definition at line 1 of file pdtrrfs.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, LIWORK, LWORK,
13 $ N, NRHS
14* ..
15* .. Array Arguments ..
16 INTEGER DESCA( * ), DESCB( * ), DESCX( * ), IWORK( * )
17 DOUBLE PRECISION A( * ), B( * ), BERR( * ), FERR( * ),
18 $ WORK( * ), X( * )
19* ..
20*
21* Purpose
22* =======
23*
24* PDTRRFS 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 PDTRTRS or some other
29* means before entering this routine. PDTRRFS 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 )**T * sub( X ) = sub( B )
102* (Conjugate transpose = Transpose)
103*
104*
105* DIAG (global input) CHARACTER*1
106* = 'N': sub( A ) is non-unit triangular;
107* = 'U': sub( A ) is unit triangular.
108*
109* N (global input) INTEGER
110* The order of the matrix sub( A ). N >= 0.
111*
112* NRHS (global input) INTEGER
113* The number of right hand sides, i.e., the number of columns
114* of the matrices sub( B ) and sub( X ). NRHS >= 0.
115*
116* A (local input) DOUBLE PRECISION pointer into the local memory
117* to an array of local dimension (LLD_A,LOCc(JA+N-1) ). This
118* array contains the local pieces of the original triangular
119* distributed matrix sub( A ).
120* If UPLO = 'U', the leading N-by-N upper triangular part of
121* sub( A ) contains the upper triangular part of the matrix,
122* and its strictly lower triangular part is not referenced.
123* If UPLO = 'L', the leading N-by-N lower triangular part of
124* sub( A ) contains the lower triangular part of the distribu-
125* ted matrix, and its strictly upper triangular part is not
126* referenced.
127* If DIAG = 'U', the diagonal elements of sub( A ) are also
128* not referenced and are assumed to be 1.
129*
130* IA (global input) INTEGER
131* The row index in the global array A indicating the first
132* row of sub( A ).
133*
134* JA (global input) INTEGER
135* The column index in the global array A indicating the
136* first column of sub( A ).
137*
138* DESCA (global and local input) INTEGER array of dimension DLEN_.
139* The array descriptor for the distributed matrix A.
140*
141* B (local input) DOUBLE PRECISION pointer into the local memory
142* to an array of local dimension (LLD_B, LOCc(JB+NRHS-1) ).
143* On entry, this array contains the the local pieces of the
144* right hand sides sub( B ).
145*
146* IB (global input) INTEGER
147* The row index in the global array B indicating the first
148* row of sub( B ).
149*
150* JB (global input) INTEGER
151* The column index in the global array B indicating the
152* first column of sub( B ).
153*
154* DESCB (global and local input) INTEGER array of dimension DLEN_.
155* The array descriptor for the distributed matrix B.
156*
157* X (local input) DOUBLE PRECISION pointer into the local memory
158* to an array of local dimension (LLD_X, LOCc(JX+NRHS-1) ).
159* On entry, this array contains the the local pieces of the
160* solution vectors sub( X ).
161*
162* IX (global input) INTEGER
163* The row index in the global array X indicating the first
164* row of sub( X ).
165*
166* JX (global input) INTEGER
167* The column index in the global array X indicating the
168* first column of sub( X ).
169*
170* DESCX (global and local input) INTEGER array of dimension DLEN_.
171* The array descriptor for the distributed matrix X.
172*
173* FERR (local output) DOUBLE PRECISION array of local dimension
174* LOCc(JB+NRHS-1). The estimated forward error bounds for
175* each solution vector of sub( X ). If XTRUE is the true
176* solution, FERR bounds the magnitude of the largest entry
177* in (sub( X ) - XTRUE) divided by the magnitude of the
178* largest entry in sub( X ). The estimate is as reliable as
179* the estimate for RCOND, and is almost always a slight
180* overestimate of the true error.
181* This array is tied to the distributed matrix X.
182*
183* BERR (local output) DOUBLE PRECISION array of local dimension
184* LOCc(JB+NRHS-1). The componentwise relative backward
185* error of each solution vector (i.e., the smallest re-
186* lative change in any entry of sub( A ) or sub( B )
187* that makes sub( X ) an exact solution).
188* This array is tied to the distributed matrix X.
189*
190* WORK (local workspace/local output) DOUBLE PRECISION array,
191* dimension (LWORK)
192* On exit, WORK(1) returns the minimal and optimal LWORK.
193*
194* LWORK (local or global input) INTEGER
195* The dimension of the array WORK.
196* LWORK is local input and must be at least
197* LWORK >= 3*LOCr( N + MOD( IA-1, MB_A ) ).
198*
199* If LWORK = -1, then LWORK is global input and a workspace
200* query is assumed; the routine only calculates the minimum
201* and optimal size for all work arrays. Each of these
202* values is returned in the first entry of the corresponding
203* work array, and no error message is issued by PXERBLA.
204*
205* IWORK (local workspace/local output) INTEGER array,
206* dimension (LIWORK)
207* On exit, IWORK(1) returns the minimal and optimal LIWORK.
208*
209* LIWORK (local or global input) INTEGER
210* The dimension of the array IWORK.
211* LIWORK is local input and must be at least
212* LIWORK >= LOCr( N + MOD( IB-1, MB_B ) ).
213*
214* If LIWORK = -1, then LIWORK is global input and a workspace
215* query is assumed; the routine only calculates the minimum
216* and optimal size for all work arrays. Each of these
217* values is returned in the first entry of the corresponding
218* work array, and no error message is issued by PXERBLA.
219*
220*
221* INFO (global output) INTEGER
222* = 0: successful exit
223* < 0: If the i-th argument is an array and the j-entry had
224* an illegal value, then INFO = -(i*100+j), if the i-th
225* argument is a scalar and had an illegal value, then
226* INFO = -i.
227*
228* Notes
229* =====
230*
231* This routine temporarily returns when N <= 1.
232*
233* The distributed submatrices sub( X ) and sub( B ) should be
234* distributed the same way on the same processes. These conditions
235* ensure that sub( X ) and sub( B ) are "perfectly" aligned.
236*
237* Moreover, this routine requires the distributed submatrices sub( A ),
238* sub( X ), and sub( B ) to be aligned on a block boundary,
239* i.e., if f(x,y) = MOD( x-1, y ):
240* f( IA, DESCA( MB_ ) ) = f( JA, DESCA( NB_ ) ) = 0,
241* f( IB, DESCB( MB_ ) ) = f( JB, DESCB( NB_ ) ) = 0, and
242* f( IX, DESCX( MB_ ) ) = f( JX, DESCX( NB_ ) ) = 0.
243*
244* =====================================================================
245*
246* .. Parameters ..
247 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
248 $ LLD_, MB_, M_, NB_, N_, RSRC_
249 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
250 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
251 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
252 DOUBLE PRECISION ZERO, ONE
253 parameter( zero = 0.0d+0, one = 1.0d+0 )
254* ..
255* .. Local Scalars ..
256 LOGICAL LQUERY, NOTRAN, NOUNIT, UPPER
257 CHARACTER TRANST
258 INTEGER IAROW, IXBCOL, IXBROW, IXCOL, IXROW, ICOFFA,
259 $ ICOFFB, ICOFFX, ICTXT, ICURCOL, IDUM, II, IIXB,
260 $ IIW, IOFFXB, IPB, IPR, IPV, IROFFA, IROFFB,
261 $ IROFFX, IW, J, JBRHS, JJ, JJFBE, JJXB, JN, JW,
262 $ K, KASE, LDXB, LIWMIN, LWMIN, MYCOL, MYRHS,
263 $ MYROW, NP, NP0, NPCOL, NPMOD, NPROW, NZ
264 DOUBLE PRECISION EPS, EST, LSTRES, S, SAFE1, SAFE2, SAFMIN
265* ..
266* .. Local Arrays ..
267 INTEGER DESCW( DLEN_ ), IDUM1( 5 ), IDUM2( 5 )
268* ..
269* .. External Functions ..
270 LOGICAL LSAME
271 INTEGER ICEIL, INDXG2P, NUMROC
272 DOUBLE PRECISION PDLAMCH
273 EXTERNAL iceil, indxg2p, lsame, numroc, pdlamch
274* ..
275* .. External Subroutines ..
276 EXTERNAL blacs_gridinfo, chk1mat, descset, dgamx2d,
277 $ dgebr2d, dgebs2d, infog2l, pchk1mat,
278 $ pchk2mat, pdatrmv, pdaxpy, pdcopy,
279 $ pdlacon, pdtrmv, pdtrsv, pxerbla
280* ..
281* .. Intrinsic Functions ..
282 INTRINSIC abs, dble, ichar, max, min, mod
283* ..
284* .. Executable Statements ..
285*
286* Get grid parameters
287*
288 ictxt = desca( ctxt_ )
289 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
290*
291* Test the input parameters.
292*
293 info = 0
294 IF( nprow.EQ.-1 ) THEN
295 info = -( 900+ctxt_ )
296 ELSE
297 CALL chk1mat( n, 4, n, 4, ia, ja, desca, 9, info )
298 CALL chk1mat( n, 4, nrhs, 5, ib, jb, descb, 13, info )
299 CALL chk1mat( n, 4, nrhs, 5, ix, jx, descx, 17, info )
300 IF( info.EQ.0 ) THEN
301 upper = lsame( uplo, 'U' )
302 notran = lsame( trans, 'N' )
303 nounit = lsame( diag, 'N' )
304 iroffa = mod( ia-1, desca( mb_ ) )
305 icoffa = mod( ja-1, desca( nb_ ) )
306 iroffb = mod( ib-1, descb( mb_ ) )
307 icoffb = mod( jb-1, descb( nb_ ) )
308 iroffx = mod( ix-1, descx( mb_ ) )
309 icoffx = mod( jx-1, descx( nb_ ) )
310 iarow = indxg2p( ia, desca( mb_ ), myrow, desca( rsrc_ ),
311 $ nprow )
312 CALL infog2l( ib, jb, descb, nprow, npcol, myrow, mycol,
313 $ iixb, jjxb, ixbrow, ixbcol )
314 ixrow = indxg2p( ix, descx( mb_ ), myrow, descx( rsrc_ ),
315 $ nprow )
316 ixcol = indxg2p( jx, descx( nb_ ), mycol, descx( csrc_ ),
317 $ npcol )
318 npmod = numroc( n+iroffa, desca( mb_ ), myrow, iarow,
319 $ nprow )
320 lwmin = 3*npmod
321 work( 1 ) = dble( lwmin )
322 liwmin = npmod
323 iwork( 1 ) = liwmin
324 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
325*
326 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
327 info = -1
328 ELSE IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND.
329 $ .NOT.lsame( trans, 'C' ) ) THEN
330 info = -2
331 ELSE IF( .NOT.nounit .AND. .NOT.lsame( diag, 'U' ) ) THEN
332 info = -3
333 ELSE IF( n.LT.0 ) THEN
334 info = -4
335 ELSE IF( nrhs.LT.0 ) THEN
336 info = -5
337 ELSE IF( iroffa.NE.0 ) THEN
338 info = -7
339 ELSE IF( icoffa.NE.0 ) THEN
340 info = -8
341 ELSE IF( desca( mb_ ).NE.desca( nb_ ) ) THEN
342 info = -( 900+nb_ )
343 ELSE IF( iroffa.NE.iroffb .OR. iarow.NE.ixbrow ) THEN
344 info = -11
345 ELSE IF( desca( mb_ ).NE.descb( mb_ ) ) THEN
346 info = -( 1300+mb_ )
347 ELSE IF( ictxt.NE.descb( ctxt_ ) ) THEN
348 info = -( 1300+ctxt_ )
349 ELSE IF( iroffx.NE.0 .OR. ixbrow.NE.ixrow ) THEN
350 info = -15
351 ELSE IF( icoffb.NE.icoffx .OR. ixbcol.NE.ixcol ) THEN
352 info = -16
353 ELSE IF( descb( mb_ ).NE.descx( mb_ ) ) THEN
354 info = -( 1700+mb_ )
355 ELSE IF( descb( nb_ ).NE.descx( nb_ ) ) THEN
356 info = -( 1700+nb_ )
357 ELSE IF( ictxt.NE.descx( ctxt_ ) ) THEN
358 info = -( 1700+ctxt_ )
359 ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
360 info = -21
361 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
362 info = -23
363 END IF
364 END IF
365*
366 IF( upper ) THEN
367 idum1( 1 ) = ichar( 'U' )
368 ELSE
369 idum1( 1 ) = ichar( 'L' )
370 END IF
371 idum2( 1 ) = 1
372 IF( notran ) THEN
373 idum1( 2 ) = ichar( 'N' )
374 ELSE IF( lsame( trans, 'T' ) ) THEN
375 idum1( 2 ) = ichar( 'T' )
376 ELSE
377 idum1( 2 ) = ichar( 'C' )
378 END IF
379 idum2( 2 ) = 2
380 IF( nounit ) THEN
381 idum1( 3 ) = ichar( 'N' )
382 ELSE
383 idum1( 3 ) = ichar( 'U' )
384 END IF
385 idum2( 3 ) = 3
386 IF( lwork.EQ.-1 ) THEN
387 idum1( 4 ) = -1
388 ELSE
389 idum1( 4 ) = 1
390 END IF
391 idum2( 4 ) = 21
392 IF( liwork.EQ.-1 ) THEN
393 idum1( 5 ) = -1
394 ELSE
395 idum1( 5 ) = 1
396 END IF
397 idum2( 5 ) = 23
398 CALL pchk1mat( n, 4, n, 4, ia, ja, desca, 9, 0, idum1, idum2,
399 $ info )
400 CALL pchk2mat( n, 4, nrhs, 5, ib, jb, descb, 13, n, 4, nrhs, 5,
401 $ ix, jx, descx, 17, 5, idum1, idum2, info )
402 END IF
403 IF( info.NE.0 ) THEN
404 CALL pxerbla( ictxt, 'PDTRRFS', -info )
405 RETURN
406 ELSE IF( lquery ) THEN
407 RETURN
408 END IF
409*
410 jjfbe = jjxb
411 myrhs = numroc( jb+nrhs-1, descb( nb_ ), mycol, descb( csrc_ ),
412 $ npcol )
413*
414* Quick return if possible
415*
416 IF( n.LE.1 .OR. nrhs.EQ.0 ) THEN
417 DO 10 jj = jjfbe, myrhs
418 ferr( jj ) = zero
419 berr( jj ) = zero
420 10 CONTINUE
421 RETURN
422 END IF
423*
424 IF( notran ) THEN
425 transt = 'T'
426 ELSE
427 transt = 'N'
428 END IF
429*
430 np0 = numroc( n+iroffb, descb( mb_ ), myrow, ixbrow, nprow )
431 CALL descset( descw, n+iroffb, 1, desca( mb_ ), 1, ixbrow, ixbcol,
432 $ ictxt, max( 1, np0 ) )
433 ipb = 1
434 ipr = ipb + np0
435 ipv = ipr + np0
436 IF( myrow.EQ.ixbrow ) THEN
437 iiw = 1 + iroffb
438 np = np0 - iroffb
439 ELSE
440 iiw = 1
441 np = np0
442 END IF
443 iw = 1 + iroffb
444 jw = 1
445 ldxb = descb( lld_ )
446 ioffxb = ( jjxb-1 )*ldxb
447*
448* NZ = maximum number of nonzero entries in each row of A, plus 1
449*
450 nz = n + 1
451 eps = pdlamch( ictxt, 'Epsilon' )
452 safmin = pdlamch( ictxt, 'Safe minimum' )
453 safe1 = nz*safmin
454 safe2 = safe1 / eps
455 jn = min( iceil( jb, descb( nb_ ) )*descb( nb_ ), jb+nrhs-1 )
456*
457* Handle first block separately
458*
459 jbrhs = jn - jb + 1
460 DO 90 k = 0, jbrhs - 1
461*
462* Compute residual R = B - op(A) * X,
463* where op(A) = A or A', depending on TRANS.
464*
465 CALL pdcopy( n, x, ix, jx+k, descx, 1, work( ipr ), iw, jw,
466 $ descw, 1 )
467 CALL pdtrmv( uplo, trans, diag, n, a, ia, ja, desca,
468 $ work( ipr ), iw, jw, descw, 1 )
469 CALL pdaxpy( n, -one, b, ib, jb+k, descb, 1, work( ipr ), iw,
470 $ jw, descw, 1 )
471*
472* Compute componentwise relative backward error from formula
473*
474* max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
475*
476* where abs(Z) is the componentwise absolute value of the matrix
477* or vector Z. If the i-th component of the denominator is less
478* than SAFE2, then SAFE1 is added to the i-th components of the
479* numerator and denominator before dividing.
480*
481 IF( mycol.EQ.ixbcol ) THEN
482 IF( np.GT.0 ) THEN
483 DO 20 ii = iixb, iixb + np - 1
484 work( iiw+ii-iixb ) = abs( b( ii+ioffxb ) )
485 20 CONTINUE
486 END IF
487 END IF
488*
489 CALL pdatrmv( uplo, trans, diag, n, one, a, ia, ja, desca, x,
490 $ ix, jx+k, descx, 1, one, work( ipb ), iw, jw,
491 $ descw, 1 )
492*
493 s = zero
494 IF( mycol.EQ.ixbcol ) THEN
495 IF( np.GT.0 ) THEN
496 DO 30 ii = iiw - 1, iiw + np - 2
497 IF( work( ipb+ii ).GT.safe2 ) THEN
498 s = max( s, abs( work( ipr+ii ) ) /
499 $ work( ipb+ii ) )
500 ELSE
501 s = max( s, ( abs( work( ipr+ii ) )+safe1 ) /
502 $ ( work( ipb+ii )+safe1 ) )
503 END IF
504 30 CONTINUE
505 END IF
506 END IF
507*
508 CALL dgamx2d( ictxt, 'All', ' ', 1, 1, s, 1, idum, idum, 1,
509 $ -1, mycol )
510 IF( mycol.EQ.ixbcol )
511 $ berr( jjfbe ) = s
512*
513* Bound error from formula
514*
515* norm(X - XTRUE) / norm(X) .le. FERR =
516* norm( abs(inv(op(A)))*
517* ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
518*
519* where
520* norm(Z) is the magnitude of the largest component of Z
521* inv(op(A)) is the inverse of op(A)
522* abs(Z) is the componentwise absolute value of the matrix or
523* vector Z
524* NZ is the maximum number of nonzeros in any row of A, plus 1
525* EPS is machine epsilon
526*
527* The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
528* is incremented by SAFE1 if the i-th component of
529* abs(op(A))*abs(X) + abs(B) is less than SAFE2.
530*
531* Use PDLACON to estimate the infinity-norm of the matrix
532* inv(op(A)) * diag(W),
533* where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
534*
535 IF( mycol.EQ.ixbcol ) THEN
536 IF( np.GT.0 ) THEN
537 DO 40 ii = iiw - 1, iiw + np - 2
538 IF( work( ipb+ii ).GT.safe2 ) THEN
539 work( ipb+ii ) = abs( work( ipr+ii ) ) +
540 $ nz*eps*work( ipb+ii )
541 ELSE
542 work( ipb+ii ) = abs( work( ipr+ii ) ) +
543 $ nz*eps*work( ipb+ii ) + safe1
544 END IF
545 40 CONTINUE
546 END IF
547 END IF
548*
549 kase = 0
550 50 CONTINUE
551 IF( mycol.EQ.ixbcol ) THEN
552 CALL dgebs2d( ictxt, 'Rowwise', ' ', np, 1, work( ipr ),
553 $ descw( lld_ ) )
554 ELSE
555 CALL dgebr2d( ictxt, 'Rowwise', ' ', np, 1, work( ipr ),
556 $ descw( lld_ ), myrow, ixbcol )
557 END IF
558 descw( csrc_ ) = mycol
559 CALL pdlacon( n, work( ipv ), iw, jw, descw, work( ipr ),
560 $ iw, jw, descw, iwork, est, kase )
561 descw( csrc_ ) = ixbcol
562*
563 IF( kase.NE.0 ) THEN
564 IF( kase.EQ.1 ) THEN
565*
566* Multiply by diag(W)*inv(op(A)').
567*
568 CALL pdtrsv( uplo, transt, diag, n, a, ia, ja, desca,
569 $ work( ipr ), iw, jw, descw, 1 )
570 IF( mycol.EQ.ixbcol ) THEN
571 IF( np.GT.0 ) THEN
572 DO 60 ii = iiw - 1, iiw + np - 2
573 work( ipr+ii ) = work( ipb+ii )*work( ipr+ii )
574 60 CONTINUE
575 END IF
576 END IF
577 ELSE
578*
579* Multiply by inv(op(A))*diag(W).
580*
581 IF( mycol.EQ.ixbcol ) THEN
582 IF( np.GT.0 ) THEN
583 DO 70 ii = iiw - 1, iiw + np - 2
584 work( ipr+ii ) = work( ipb+ii )*work( ipr+ii )
585 70 CONTINUE
586 END IF
587 END IF
588 CALL pdtrsv( uplo, trans, diag, n, a, ia, ja, desca,
589 $ work( ipr ), iw, jw, descw, 1 )
590 END IF
591 GO TO 50
592 END IF
593*
594* Normalize error.
595*
596 lstres = zero
597 IF( mycol.EQ.ixbcol ) THEN
598 IF( np.GT.0 ) THEN
599 DO 80 ii = iixb, iixb + np - 1
600 lstres = max( lstres, abs( x( ioffxb+ii ) ) )
601 80 CONTINUE
602 END IF
603 CALL dgamx2d( ictxt, 'Column', ' ', 1, 1, lstres, 1, idum,
604 $ idum, 1, -1, mycol )
605 IF( lstres.NE.zero )
606 $ ferr( jjfbe ) = est / lstres
607*
608 jjxb = jjxb + 1
609 jjfbe = jjfbe + 1
610 ioffxb = ioffxb + ldxb
611*
612 END IF
613*
614 90 CONTINUE
615*
616 icurcol = mod( ixbcol+1, npcol )
617*
618* Do for each right hand side
619*
620 DO 180 j = jn + 1, jb + nrhs - 1, descb( nb_ )
621 jbrhs = min( jb+nrhs-j, descb( nb_ ) )
622 descw( csrc_ ) = icurcol
623*
624 DO 170 k = 0, jbrhs - 1
625*
626* Compute residual R = B - op(A) * X,
627* where op(A) = A or A', depending on TRANS.
628*
629 CALL pdcopy( n, x, ix, j+k, descx, 1, work( ipr ), iw, jw,
630 $ descw, 1 )
631 CALL pdtrmv( uplo, trans, diag, n, a, ia, ja, desca,
632 $ work( ipr ), iw, jw, descw, 1 )
633 CALL pdaxpy( n, -one, b, ib, j+k, descb, 1, work( ipr ),
634 $ iw, jw, descw, 1 )
635*
636* Compute componentwise relative backward error from formula
637*
638* max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
639*
640* where abs(Z) is the componentwise absolute value of the
641* matrix or vector Z. If the i-th component of the
642* denominator is less than SAFE2, then SAFE1 is added to the
643* i-th components of the numerator and denominator before
644* dividing.
645*
646 IF( mycol.EQ.ixbcol ) THEN
647 IF( np.GT.0 ) THEN
648 DO 100 ii = iixb, iixb + np - 1
649 work( iiw+ii-iixb ) = abs( b( ii+ioffxb ) )
650 100 CONTINUE
651 END IF
652 END IF
653*
654 CALL pdatrmv( uplo, trans, diag, n, one, a, ia, ja, desca,
655 $ x, ix, j+k, descx, 1, one, work( ipb ), iw,
656 $ jw, descw, 1 )
657*
658 s = zero
659 IF( mycol.EQ.ixbcol ) THEN
660 IF( np.GT.0 ) THEN
661 DO 110 ii = iiw - 1, iiw + np - 2
662 IF( work( ipb+ii ).GT.safe2 ) THEN
663 s = max( s, abs( work( ipr+ii ) ) /
664 $ work( ipb+ii ) )
665 ELSE
666 s = max( s, ( abs( work( ipr+ii ) )+safe1 ) /
667 $ ( work( ipb+ii )+safe1 ) )
668 END IF
669 110 CONTINUE
670 END IF
671 END IF
672*
673 CALL dgamx2d( ictxt, 'All', ' ', 1, 1, s, 1, idum, idum, 1,
674 $ -1, mycol )
675 IF( mycol.EQ.ixbcol )
676 $ berr( jjfbe ) = s
677*
678* Bound error from formula
679*
680* norm(X - XTRUE) / norm(X) .le. FERR =
681* norm( abs(inv(op(A)))*
682* ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))/norm(X)
683*
684* where
685* norm(Z) is the magnitude of the largest component of Z
686* inv(op(A)) is the inverse of op(A)
687* abs(Z) is the componentwise absolute value of the matrix
688* or vector Z
689* NZ is the maximum number of nonzeros in any row of A,
690* plus 1
691* EPS is machine epsilon
692*
693* The i-th component of
694* abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
695* is incremented by SAFE1 if the i-th component of
696* abs(op(A))*abs(X) + abs(B) is less than SAFE2.
697*
698* Use PDLACON to estimate the infinity-norm of the matrix
699* inv(op(A)) * diag(W),
700* where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
701*
702 IF( mycol.EQ.ixbcol ) THEN
703 IF( np.GT.0 ) THEN
704 DO 120 ii = iiw - 1, iiw + np - 2
705 IF( work( ipb+ii ).GT.safe2 ) THEN
706 work( ipb+ii ) = abs( work( ipr+ii ) ) +
707 $ nz*eps*work( ipb+ii )
708 ELSE
709 work( ipb+ii ) = abs( work( ipr+ii ) ) +
710 $ nz*eps*work( ipb+ii ) + safe1
711 END IF
712 120 CONTINUE
713 END IF
714 END IF
715*
716 kase = 0
717 130 CONTINUE
718 IF( mycol.EQ.ixbcol ) THEN
719 CALL dgebs2d( ictxt, 'Rowwise', ' ', np, 1, work( ipr ),
720 $ descw( lld_ ) )
721 ELSE
722 CALL dgebr2d( ictxt, 'Rowwise', ' ', np, 1, work( ipr ),
723 $ descw( lld_ ), myrow, ixbcol )
724 END IF
725 descw( csrc_ ) = mycol
726 CALL pdlacon( n, work( ipv ), iw, jw, descw, work( ipr ),
727 $ iw, jw, descw, iwork, est, kase )
728 descw( csrc_ ) = ixbcol
729*
730 IF( kase.NE.0 ) THEN
731 IF( kase.EQ.1 ) THEN
732*
733* Multiply by diag(W)*inv(op(A)').
734*
735 CALL pdtrsv( uplo, transt, diag, n, a, ia, ja, desca,
736 $ work( ipr ), iw, jw, descw, 1 )
737 IF( mycol.EQ.ixbcol ) THEN
738 IF( np.GT.0 ) THEN
739 DO 140 ii = iiw - 1, iiw + np - 2
740 work( ipr+ii ) = work( ipb+ii )*
741 $ work( ipr+ii )
742 140 CONTINUE
743 END IF
744 END IF
745 ELSE
746*
747* Multiply by inv(op(A))*diag(W).
748*
749 IF( mycol.EQ.ixbcol ) THEN
750 IF( np.GT.0 ) THEN
751 DO 150 ii = iiw - 1, iiw + np - 2
752 work( ipr+ii ) = work( ipb+ii )*
753 $ work( ipr+ii )
754 150 CONTINUE
755 END IF
756 END IF
757 CALL pdtrsv( uplo, trans, diag, n, a, ia, ja, desca,
758 $ work( ipr ), iw, jw, descw, 1 )
759 END IF
760 GO TO 130
761 END IF
762*
763* Normalize error.
764*
765 lstres = zero
766 IF( mycol.EQ.ixbcol ) THEN
767 IF( np.GT.0 ) THEN
768 DO 160 ii = iixb, iixb + np - 1
769 lstres = max( lstres, abs( x( ioffxb+ii ) ) )
770 160 CONTINUE
771 END IF
772 CALL dgamx2d( ictxt, 'Column', ' ', 1, 1, lstres, 1,
773 $ idum, idum, 1, -1, mycol )
774 IF( lstres.NE.zero )
775 $ ferr( jjfbe ) = est / lstres
776*
777 jjxb = jjxb + 1
778 jjfbe = jjfbe + 1
779 ioffxb = ioffxb + ldxb
780*
781 END IF
782*
783 170 CONTINUE
784*
785 icurcol = mod( icurcol+1, npcol )
786*
787 180 CONTINUE
788*
789 work( 1 ) = dble( lwmin )
790 iwork( 1 ) = liwmin
791*
792 RETURN
793*
794* End of PDTRRFS
795*
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
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
pdlamch
double precision function pdlamch(ictxt, cmach)
Definition pdblastst.f:6769
pdlacon
subroutine pdlacon(n, v, iv, jv, descv, x, ix, jx, descx, isgn, est, kase)
Definition pdlacon.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: