ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pslaschk.f
Go to the documentation of this file.
1  SUBROUTINE pslaschk( SYMM, DIAG, N, NRHS, X, IX, JX, DESCX,
2  $ IASEED, IA, JA, DESCA, IBSEED, ANORM, RESID,
3  $ WORK )
4 *
5 * -- ScaLAPACK auxiliary 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, SYMM
12  INTEGER IA, IASEED, IBSEED, IX, JA, JX, N, NRHS
13  REAL ANORM, RESID
14 * ..
15 * .. Array Arguments ..
16  INTEGER DESCA( * ), DESCX( * )
17  REAL WORK( * ), X( * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * PSLASCHK computes the residual
24 * || sub( A )*sub( X ) - B || / (|| sub( A ) ||*|| sub( X ) ||*eps*N)
25 * to check the accuracy of the factorization and solve steps in the
26 * LU and Cholesky decompositions, where sub( A ) denotes
27 * A(IA:IA+N-1,JA,JA+N-1), sub( X ) denotes X(IX:IX+N-1, JX:JX+NRHS-1).
28 *
29 * Notes
30 * =====
31 *
32 * Each global data object is described by an associated description
33 * vector. This vector stores the information required to establish
34 * the mapping between an object element and its corresponding process
35 * and memory location.
36 *
37 * Let A be a generic term for any 2D block cyclicly distributed array.
38 * Such a global array has an associated description vector DESCA.
39 * In the following comments, the character _ should be read as
40 * "of the global array".
41 *
42 * NOTATION STORED IN EXPLANATION
43 * --------------- -------------- --------------------------------------
44 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
45 * DTYPE_A = 1.
46 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
47 * the BLACS process grid A is distribu-
48 * ted over. The context itself is glo-
49 * bal, but the handle (the integer
50 * value) may vary.
51 * M_A (global) DESCA( M_ ) The number of rows in the global
52 * array A.
53 * N_A (global) DESCA( N_ ) The number of columns in the global
54 * array A.
55 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
56 * the rows of the array.
57 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
58 * the columns of the array.
59 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
60 * row of the array A is distributed.
61 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the
62 * first column of the array A is
63 * distributed.
64 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
65 * array. LLD_A >= MAX(1,LOCr(M_A)).
66 *
67 * Let K be the number of rows or columns of a distributed matrix,
68 * and assume that its process grid has dimension p x q.
69 * LOCr( K ) denotes the number of elements of K that a process
70 * would receive if K were distributed over the p processes of its
71 * process column.
72 * Similarly, LOCc( K ) denotes the number of elements of K that a
73 * process would receive if K were distributed over the q processes of
74 * its process row.
75 * The values of LOCr() and LOCc() may be determined via a call to the
76 * ScaLAPACK tool function, NUMROC:
77 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
78 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
79 * An upper bound for these quantities may be computed by:
80 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
81 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
82 *
83 * Arguments
84 * =========
85 *
86 * SYMM (global input) CHARACTER
87 * if SYMM = 'S', sub( A ) is a symmetric distributed matrix,
88 * otherwise sub( A ) is a general distributed matrix.
89 *
90 * DIAG (global input) CHARACTER
91 * If DIAG = 'D', sub( A ) is diagonally dominant.
92 *
93 * N (global input) INTEGER
94 * The number of columns to be operated on, i.e. the number of
95 * columns of the distributed submatrix sub( A ). N >= 0.
96 *
97 * NRHS (global input) INTEGER
98 * The number of right-hand-sides, i.e the number of columns
99 * of the distributed matrix sub( X ). NRHS >= 0.
100 *
101 * X (local input) REAL pointer into the local memory
102 * to an array of dimension (LLD_X,LOCc(JX+NRHS-1). This array
103 * contains the local pieces of the answer vector(s) sub( X ) of
104 * sub( A ) sub( X ) - B, split up over a column of processes.
105 *
106 * IX (global input) INTEGER
107 * The row index in the global array X indicating the first
108 * row of sub( X ).
109 *
110 * JX (global input) INTEGER
111 * The column index in the global array X indicating the
112 * first column of sub( X ).
113 *
114 * DESCX (global and local input) INTEGER array of dimension DLEN_.
115 * The array descriptor for the distributed matrix X.
116 *
117 * IASEED (global input) INTEGER
118 * The seed number to generate the original matrix Ao.
119 *
120 * IA (global input) INTEGER
121 * The row index in the global array A indicating the first
122 * row of sub( A ).
123 *
124 * JA (global input) INTEGER
125 * The column index in the global array A indicating the
126 * first column of sub( A ).
127 *
128 * DESCA (global and local input) INTEGER array of dimension DLEN_.
129 * The array descriptor for the distributed matrix A.
130 *
131 * IBSEED (global input) INTEGER
132 * The seed number to generate the original matrix B.
133 *
134 * ANORM (global input) REAL
135 * The 1-norm or infinity norm of the distributed matrix
136 * sub( A ).
137 *
138 * RESID (global output) REAL
139 * The residual error:
140 * ||sub( A )*sub( X )-B|| / (||sub( A )||*||sub( X )||*eps*N).
141 *
142 * WORK (local workspace) REAL array, dimension (LWORK)
143 * LWORK >= MAX(1,Np)*NB_X + Nq*NB_X + MAX( MAX(NQ*MB_A,2*NB_X),
144 * NB_X * NUMROC( NUMROC(N,MB_X,0,0,NPCOL), MB_X, 0, 0, LCMQ ) )
145 *
146 * =====================================================================
147 *
148 * .. Parameters ..
149  INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
150  $ LLD_, MB_, M_, NB_, N_, RSRC_
151  parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
152  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
153  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
154  REAL ZERO, ONE
155  PARAMETER ( ONE = 1.0e+0, zero = 0.0e+0 )
156 * ..
157 * .. Local Scalars ..
158  INTEGER IACOL, IAROW, IB, ICOFF, ICTXT, ICURCOL, IDUMM,
159  $ II, IIA, IIX, IOFFX, IPA, IPB, IPW, IPX, IROFF,
160  $ ixcol, ixrow, j, jbrhs, jj, jja, jjx, ldx,
161  $ mycol, myrow, np, npcol, nprow, nq
162  REAL BETA, DIVISOR, EPS, RESID1
163 * ..
164 * .. External Subroutines ..
165  EXTERNAL blacs_gridinfo, pbstran, psmatgen,
166  $ sgamx2d, sgebr2d, sgebs2d, sgemm,
167  $ sgerv2d, sgesd2d, sgsum2d, slaset
168 * ..
169 * .. External Functions ..
170  INTEGER ISAMAX, NUMROC
171  REAL PSLAMCH
172  EXTERNAL isamax, numroc, pslamch
173 * ..
174 * .. Intrinsic Functions ..
175  INTRINSIC abs, max, min, mod, real
176 * ..
177 * .. Executable Statements ..
178 *
179 * Get needed initial parameters
180 *
181  ictxt = desca( ctxt_ )
182  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
183 *
184  eps = pslamch( ictxt, 'eps' )
185  resid = 0.0e+0
186  divisor = anorm * eps * real( n )
187 *
188  CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, iia, jja,
189  $ iarow, iacol )
190  CALL infog2l( ix, jx, descx, nprow, npcol, myrow, mycol, iix, jjx,
191  $ ixrow, ixcol )
192  iroff = mod( ia-1, desca( mb_ ) )
193  icoff = mod( ja-1, desca( nb_ ) )
194  np = numroc( n+iroff, desca( mb_ ), myrow, iarow, nprow )
195  nq = numroc( n+icoff, desca( nb_ ), mycol, iacol, npcol )
196 *
197  ldx = max( 1, np )
198  ipb = 1
199  ipx = ipb + np * descx( nb_ )
200  ipa = ipx + nq * descx( nb_ )
201 *
202  IF( myrow.EQ.iarow )
203  $ np = np - iroff
204  IF( mycol.EQ.iacol )
205  $ nq = nq - icoff
206 *
207  icurcol = ixcol
208 *
209 * Loop over the rhs
210 *
211  DO 40 j = 1, nrhs, descx( nb_ )
212  jbrhs = min( descx( nb_ ), nrhs-j+1 )
213 *
214 * Transpose x from ICURCOL to all rows
215 *
216  ioffx = iix + ( jjx - 1 ) * descx( lld_ )
217  CALL pbstran( ictxt, 'Column', 'Transpose', n, jbrhs,
218  $ descx( mb_ ), x( ioffx ), descx( lld_ ), zero,
219  $ work( ipx ), jbrhs, ixrow, icurcol, -1, iacol,
220  $ work( ipa ) )
221 *
222 * Regenerate B in IXCOL
223 *
224  IF( mycol.EQ.icurcol ) THEN
225  CALL psmatgen( ictxt, 'N', 'N', descx( m_ ), descx( n_ ),
226  $ descx( mb_ ), descx( nb_ ), work( ipb ), ldx,
227  $ ixrow, ixcol, ibseed, iix-1, np, jjx-1,
228  $ jbrhs, myrow, mycol, nprow, npcol )
229  beta = one
230  ELSE
231  beta = zero
232  END IF
233 *
234  IF( nq.GT.0 ) THEN
235  DO 10 ii = iia, iia+np-1, desca( mb_ )
236  ib = min( desca( mb_ ), iia+np-ii )
237 *
238 * Regenerate ib rows of the matrix A(IA:IA+N-1,JA:JA+N-1).
239 *
240  CALL psmatgen( ictxt, symm, diag, desca( m_ ),
241  $ desca( n_ ), desca( mb_ ), desca( nb_ ),
242  $ work( ipa ), ib, desca( rsrc_ ),
243  $ desca( csrc_ ), iaseed, ii-1, ib,
244  $ jja-1, nq, myrow, mycol, nprow, npcol )
245 *
246 * Compute B <= B - A * X.
247 *
248  CALL sgemm( 'No transpose', 'Transpose', ib, jbrhs, nq,
249  $ -one, work( ipa ), ib, work( ipx ), jbrhs,
250  $ beta, work( ipb+ii-iia ), ldx )
251 *
252  10 CONTINUE
253 *
254  ELSE IF( mycol.NE.icurcol ) THEN
255 *
256  CALL slaset( 'All', np, jbrhs, zero, zero, work( ipb ),
257  $ ldx )
258 *
259  END IF
260 *
261 * Add B rowwise to ICURCOL
262 *
263  CALL sgsum2d( ictxt, 'Row', ' ', np, jbrhs, work( ipb ), ldx,
264  $ myrow, icurcol )
265 *
266  IF( mycol.EQ.icurcol ) THEN
267 *
268 * Figure || A * X - B || & || X ||
269 *
270  ipw = ipa + jbrhs
271  DO 20 jj = 0, jbrhs - 1
272  IF( np.GT.0 ) THEN
273  ii = isamax( np, work( ipb+jj*ldx ), 1 )
274  work( ipa+jj ) = abs( work( ipb+ii-1+jj*ldx ) )
275  work( ipw+jj ) = abs( x( ioffx + isamax( np,
276  $ x( ioffx + jj*descx( lld_ ) ), 1 )-1+jj*
277  $ descx( lld_ ) ) )
278  ELSE
279  work( ipa+jj ) = zero
280  work( ipw+jj ) = zero
281  END IF
282  20 CONTINUE
283 *
284 * After SGAMX2D computation,
285 * WORK(IPB) has the maximum of || Ax - b ||, and
286 * WORK(IPX) has the maximum of || X ||.
287 *
288  CALL sgamx2d( ictxt, 'Column', ' ', 1, 2*jbrhs,
289  $ work( ipa ), 1, idumm, idumm, -1, 0, icurcol )
290 *
291 * Calculate residual = ||Ax-b|| / (||x||*||A||*eps*N)
292 *
293  IF( myrow.EQ.0 ) THEN
294  DO 30 jj = 0, jbrhs - 1
295  resid1 = work( ipa+jj ) / ( work( ipw+jj )*divisor )
296  IF( resid.LT.resid1 )
297  $ resid = resid1
298  30 CONTINUE
299  IF( mycol.NE.0 )
300  $ CALL sgesd2d( ictxt, 1, 1, resid, 1, 0, 0 )
301  END IF
302 *
303  ELSE IF( myrow.EQ.0 .AND. mycol.EQ.0 ) THEN
304 *
305  CALL sgerv2d( ictxt, 1, 1, resid1, 1, 0, icurcol )
306  IF( resid.LT.resid1 )
307  $ resid = resid1
308 *
309  END IF
310 *
311  IF( mycol.EQ.icurcol )
312  $ jjx = jjx + jbrhs
313  icurcol = mod( icurcol+1, npcol )
314 *
315  40 CONTINUE
316 *
317  IF( myrow.EQ.0 .AND. mycol.EQ.0 ) THEN
318  CALL sgebs2d( ictxt, 'All', ' ', 1, 1, resid, 1 )
319  ELSE
320  CALL sgebr2d( ictxt, 'All', ' ', 1, 1, resid, 1, 0, 0 )
321  END IF
322 *
323  RETURN
324 *
325 * End of PSLASCHK
326 *
327  END
max
#define max(A, B)
Definition: pcgemr.c:180
infog2l
subroutine infog2l(GRINDX, GCINDX, DESC, NPROW, NPCOL, MYROW, MYCOL, LRINDX, LCINDX, RSRC, CSRC)
Definition: infog2l.f:3
psmatgen
subroutine psmatgen(ICTXT, AFORM, DIAG, M, N, MB, NB, A, LDA, IAROW, IACOL, ISEED, IROFF, IRNUM, ICOFF, ICNUM, MYROW, MYCOL, NPROW, NPCOL)
Definition: psmatgen.f:4
pbstran
subroutine pbstran(ICONTXT, ADIST, TRANS, M, N, NB, A, LDA, BETA, C, LDC, IAROW, IACOL, ICROW, ICCOL, WORK)
Definition: pbstran.f:3
pslaschk
subroutine pslaschk(SYMM, DIAG, N, NRHS, X, IX, JX, DESCX, IASEED, IA, JA, DESCA, IBSEED, ANORM, RESID, WORK)
Definition: pslaschk.f:4
min
#define min(A, B)
Definition: pcgemr.c:181