1      SUBROUTINE PDPBLASCHK( SYMM, UPLO, N, BWL, BWU, NRHS, X, IX, JX,
2     $                       DESCX, IASEED, A, IA, JA, DESCA, IBSEED,
3     $                       ANORM, RESID, WORK, WORKSIZ )
4*
5*
6*  -- ScaLAPACK routine (version 1.7) --
7*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
8*     and University of California, Berkeley.
9*     November 15, 1997
10*
11*     .. Scalar Arguments ..
12      CHARACTER          SYMM, UPLO
13      INTEGER            BWL, BWU, IA, IASEED, IBSEED,
14     $                   IX, JA, JX, N, NRHS, WORKSIZ
15      DOUBLE PRECISION   ANORM, RESID
16*     ..
17*     .. Array Arguments ..
18      INTEGER            DESCA( * ), DESCX( * )
19      DOUBLE PRECISION   A( * ), WORK( * ), X( * )
20*     .. External Functions ..
21      LOGICAL            LSAME
22*     ..
23*
24*  Purpose
25*  =======
26*
27*  PDPBLASCHK computes the residual
28*  || sub( A )*sub( X ) - B || / (|| sub( A ) ||*|| sub( X ) ||*eps*N)
29*  to check the accuracy of the factorization and solve steps in the
30*  LU and Cholesky decompositions, where sub( A ) denotes
31*  A(IA:IA+N-1,JA,JA+N-1), sub( X ) denotes X(IX:IX+N-1, JX:JX+NRHS-1).
32*
33*  Notes
34*  =====
35*
36*  Each global data object is described by an associated description
37*  vector.  This vector stores the information required to establish
38*  the mapping between an object element and its corresponding process
39*  and memory location.
40*
41*  Let A be a generic term for any 2D block cyclicly distributed array.
42*  Such a global array has an associated description vector DESCA.
43*  In the following comments, the character _ should be read as
44*  "of the global array".
45*
46*  NOTATION        STORED IN      EXPLANATION
47*  --------------- -------------- --------------------------------------
48*  DTYPE_A(global) DESCA( DTYPE_ )The descriptor type.  In this case,
49*                                 DTYPE_A = 1.
50*  CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
51*                                 the BLACS process grid A is distribu-
52*                                 ted over. The context itself is glo-
53*                                 bal, but the handle (the integer
54*                                 value) may vary.
55*  M_A    (global) DESCA( M_ )    The number of rows in the global
56*                                 array A.
57*  N_A    (global) DESCA( N_ )    The number of columns in the global
58*                                 array A.
59*  MB_A   (global) DESCA( MB_ )   The blocking factor used to distribute
60*                                 the rows of the array.
61*  NB_A   (global) DESCA( NB_ )   The blocking factor used to distribute
62*                                 the columns of the array.
63*  RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
64*                                 row of the array A is distributed.
65*  CSRC_A (global) DESCA( CSRC_ ) The process column over which the
66*                                 first column of the array A is
67*                                 distributed.
68*  LLD_A  (local)  DESCA( LLD_ )  The leading dimension of the local
69*                                 array.  LLD_A >= MAX(1,LOCr(M_A)).
70*
71*  Let K be the number of rows or columns of a distributed matrix,
72*  and assume that its process grid has dimension p x q.
73*  LOCr( K ) denotes the number of elements of K that a process
74*  would receive if K were distributed over the p processes of its
75*  process column.
76*  Similarly, LOCc( K ) denotes the number of elements of K that a
77*  process would receive if K were distributed over the q processes of
78*  its process row.
79*  The values of LOCr() and LOCc() may be determined via a call to the
80*  ScaLAPACK tool function, NUMROC:
81*          LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
82*          LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
83*  An upper bound for these quantities may be computed by:
84*          LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
85*          LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
86*
87*  Arguments
88*  =========
89*
90*  SYMM      (global input) CHARACTER
91*          if SYMM = 'S', sub( A ) is a symmetric distributed band
92*          matrix, otherwise sub( A ) is a general distributed matrix.
93*
94*  UPLO    (global input) CHARACTER
95*          if SYMM = 'S', then
96*            if UPLO = 'L', the lower half of the matrix is stored
97*            if UPLO = 'U', the upper half of the matrix is stored
98*          if SYMM != 'S' or 'H', then
99*            if UPLO = 'D', the matrix is stable during factorization
100*                           without interchanges
101*            if UPLO != 'D', the matrix is general
102*
103*  N       (global input) INTEGER
104*          The number of columns to be operated on, i.e. the number of
105*          columns of the distributed submatrix sub( A ). N >= 0.
106*
107*  NRHS    (global input) INTEGER
108*          The number of right-hand-sides, i.e the number of columns
109*          of the distributed matrix sub( X ). NRHS >= 1.
110*
111*  X       (local input) DOUBLE PRECISION pointer into the local memory
112*          to an array of dimension (LLD_X,LOCq(JX+NRHS-1). This array
113*          contains the local pieces of the answer vector(s) sub( X ) of
114*          sub( A ) sub( X ) - B, split up over a column of processes.
115*
116*  IX      (global input) INTEGER
117*          The row index in the global array X indicating the first
118*          row of sub( X ).
119*
120*  DESCX   (global and local input) INTEGER array of dimension DLEN_.
121*          The array descriptor for the distributed matrix X.
122*
123*  IASEED  (global input) INTEGER
124*          The seed number to generate the original matrix Ao.
125*
126*  JA      (global input) INTEGER
127*          The column index in the global array A indicating the
128*          first column of sub( A ).
129*
130*  DESCA   (global and local input) INTEGER array of dimension DLEN_.
131*          The array descriptor for the distributed matrix A.
132*
133*  IBSEED  (global input) INTEGER
134*          The seed number to generate the original matrix B.
135*
136*  ANORM   (global input) DOUBLE PRECISION
137*          The 1-norm or infinity norm of the distributed matrix
138*          sub( A ).
139*
140*  RESID   (global output) DOUBLE PRECISION
141*          The residual error:
142*          ||sub( A )*sub( X )-B|| / (||sub( A )||*||sub( X )||*eps*N).
143*
144*  WORK    (local workspace) DOUBLE PRECISION array, dimension (LWORK)
145*          IF SYMM='S'
146*          LWORK >= max(5,max(bw*(bw+2),NB))+2*NB
147*          IF SYMM!='S' or 'H'
148*          LWORK >= max(5,max(bw*(bw+2),NB))+2*NB
149*
150*  WORKSIZ (local input) size of WORK.
151*
152*  =====================================================================
153*
154*  Code Developer: Andrew J. Cleary, University of Tennessee.
155*    Current address: Lawrence Livermore National Labs.
156*  This version released: August, 2001.
157*
158*  =====================================================================
159*
160*     .. Parameters ..
161      DOUBLE PRECISION   ZERO, ONE
162      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
163      INTEGER            BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
164     $                   LLD_, MB_, M_, NB_, N_, RSRC_
165      PARAMETER          ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1,
166     $                     CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6,
167     $                     RSRC_ = 7, CSRC_ = 8, LLD_ = 9 )
168      INTEGER            INT_ONE
169      PARAMETER          ( INT_ONE = 1 )
170*     ..
171*     .. Local Scalars ..
172      INTEGER            IACOL, IAROW, ICTXT,
173     $                   IIA, IIX, IPB, IPW,
174     $                   IXCOL, IXROW, J, JJA, JJX, LDA,
175     $                   MYCOL, MYROW, NB, NP, NPCOL, NPROW, NQ
176      INTEGER            BW, INFO, IPPRODUCT, WORK_MIN
177      DOUBLE PRECISION   DIVISOR, EPS, RESID1, NORMX
178*     ..
179*     .. Local Arrays ..
180*     ..
181*     .. External Subroutines ..
182      EXTERNAL           BLACS_GRIDINFO, DGAMX2D, DGEBR2D,
183     $                   DGEBS2D, DGEMM, DGERV2D, DGESD2D,
184     $                   DGSUM2D, DLASET, PBDTRAN, PDMATGEN
185*     ..
186*     .. External Functions ..
187      INTEGER            IDAMAX, NUMROC
188      DOUBLE PRECISION   PDLAMCH
189      EXTERNAL           IDAMAX, NUMROC, PDLAMCH
190*     ..
191*     .. Intrinsic Functions ..
192      INTRINSIC          ABS, DBLE, MAX, MIN, MOD
193*     ..
194*     .. Executable Statements ..
195*
196*     Get needed initial parameters
197*
198      ICTXT = DESCA( CTXT_ )
199      NB = DESCA( NB_ )
200*
201      IF( LSAME( SYMM, 'S' ) ) THEN
202         BW = BWL
203         WORK_MIN = MAX(5,MAX(BW*(BW+2),NB))+2*NB
204      ELSE
205         BW = MAX(BWL, BWU)
206         WORK_MIN = MAX(5,MAX(BW*(BW+2),NB))+2*NB
207      ENDIF
208*
209      IF ( WORKSIZ .LT. WORK_MIN ) THEN
210       CALL PXERBLA( ICTXT, 'PDBLASCHK', -18 )
211       RETURN
212      END IF
213*
214      CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
215*
216      EPS = PDLAMCH( ICTXT, 'eps' )
217      RESID = 0.0D+0
218      DIVISOR = ANORM * EPS * DBLE( N )
219*
220      CALL INFOG2L( IA, JA, DESCA, NPROW, NPCOL, MYROW, MYCOL, IIA, JJA,
221     $              IAROW, IACOL )
222      CALL INFOG2L( IX, JX, DESCX, NPROW, NPCOL, MYROW, MYCOL, IIX, JJX,
223     $              IXROW, IXCOL )
224      NP = NUMROC( (BW+1), DESCA( MB_ ), MYROW, 0, NPROW )
225      NQ = NUMROC( N, DESCA( NB_ ), MYCOL, 0, NPCOL )
226*
227      IPB = 1
228      IPPRODUCT = 1 + DESCA( NB_ )
229      IPW = 1 + 2*DESCA( NB_ )
230*
231      LDA = DESCA( LLD_ )
232*
233*     Regenerate A
234*
235               IF( LSAME( SYMM, 'S' )) THEN
236                 CALL PDBMATGEN( ICTXT, UPLO, 'D', BW, BW, N, BW+1,
237     $                           DESCA( NB_ ), A, DESCA( LLD_ ), 0, 0,
238     $                           IASEED, MYROW, MYCOL, NPROW, NPCOL )
239               ELSE
240*
241                 CALL PDBMATGEN( ICTXT, 'N', UPLO, BWL, BWU, N,
242     $                           DESCA( MB_ ), DESCA( NB_ ), A,
243     $                           DESCA( LLD_ ), 0, 0, IASEED, MYROW,
244     $                           MYCOL, NPROW, NPCOL )
245               ENDIF
246*
247*     Loop over the rhs
248*
249      RESID = 0.0
250*
251      DO 40 J = 1, NRHS
252*
253*           Multiply A * current column of X
254*
255*
256            CALL PDPBDCMV( BW+1, BW, UPLO, N, A, 1, DESCA,
257     $          1, X( 1 + (J-1)*DESCX( LLD_ )), 1, DESCX,
258     $          WORK( IPPRODUCT ), WORK( IPW ), (BW+2)*BW, INFO )
259*
260*
261*           Regenerate column of B
262*
263            CALL PDMATGEN( DESCX( CTXT_ ), 'No', 'No', DESCX( M_ ),
264     $                     DESCX( N_ ), DESCX( MB_ ), DESCX( NB_ ),
265     $                     WORK( IPB ), DESCX( LLD_ ), DESCX( RSRC_ ),
266     $                     DESCX( CSRC_ ), IBSEED, 0, NQ, J-1, 1, MYCOL,
267     $                     MYROW, NPCOL, NPROW )
268*
269*           Figure || A * X - B || & || X ||
270*
271            CALL PDAXPY( N, -ONE, WORK( IPPRODUCT ), 1, 1, DESCX, 1,
272     $                   WORK( IPB ), 1, 1, DESCX, 1 )
273*
274            CALL PDNRM2( N, NORMX,
275     $           X, 1, J, DESCX, 1 )
276*
277            CALL PDNRM2( N, RESID1,
278     $           WORK( IPB ), 1, 1, DESCX, 1 )
279*
280*
281*           Calculate residual = ||Ax-b|| / (||x||*||A||*eps*N)
282*
283            RESID1 = RESID1 / ( NORMX*DIVISOR )
284*
285            RESID = MAX( RESID, RESID1 )
286*
287   40 CONTINUE
288*
289      RETURN
290*
291*     End of PDBLASCHK
292*
293      END
294