ScaLAPACK 2.1  2.1 ScaLAPACK: Scalable Linear Algebra PACKage
pcqrt16.f
Go to the documentation of this file.
1  SUBROUTINE pcqrt16( TRANS, M, N, NRHS, A, IA, JA, DESCA, X, IX,
2  \$ JX, DESCX, B, IB, JB, DESCB, RWORK, RESID )
3 *
4 * -- ScaLAPACK routine (version 1.7) --
5 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
6 * and University of California, Berkeley.
7 * May 1, 1997
8 *
9 * .. Scalar Arguments ..
10  CHARACTER TRANS
11  INTEGER IA, IB, IX, JA, JB, JX, M, N, NRHS
12  REAL RESID
13 * ..
14 * .. Array Arguments ..
15  INTEGER DESCA( * ), DESCB( * ), DESCX( * )
16  REAL RWORK( * )
17  COMPLEX A( * ), B( * ), X( * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * PCQRT16 computes the residual for a solution of a system of linear
24 * equations sub( A )*sub( X ) = B or sub( A' )*sub( X ) = B:
25 * RESID = norm(B - sub( A )*sub( X ) ) /
26 * ( max(m,n) * norm(sub( A ) ) * norm(sub( X ) ) * EPS ),
27 * where EPS is the machine epsilon, sub( A ) denotes
28 * A(IA:IA+N-1,JA,JA+N-1), and sub( X ) denotes
29 * X(IX:IX+N-1, JX:JX+NRHS-1).
30 *
31 * Notes
32 * =====
33 *
34 * Each global data object is described by an associated description
35 * vector. This vector stores the information required to establish
36 * the mapping between an object element and its corresponding process
37 * and memory location.
38 *
39 * Let A be a generic term for any 2D block cyclicly distributed array.
40 * Such a global array has an associated description vector DESCA.
41 * In the following comments, the character _ should be read as
42 * "of the global array".
43 *
44 * NOTATION STORED IN EXPLANATION
45 * --------------- -------------- --------------------------------------
46 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
47 * DTYPE_A = 1.
48 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
49 * the BLACS process grid A is distribu-
50 * ted over. The context itself is glo-
51 * bal, but the handle (the integer
52 * value) may vary.
53 * M_A (global) DESCA( M_ ) The number of rows in the global
54 * array A.
55 * N_A (global) DESCA( N_ ) The number of columns in the global
56 * array A.
57 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
58 * the rows of the array.
59 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
60 * the columns of the array.
61 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
62 * row of the array A is distributed.
63 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the
64 * first column of the array A is
65 * distributed.
66 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
67 * array. LLD_A >= MAX(1,LOCr(M_A)).
68 *
69 * Let K be the number of rows or columns of a distributed matrix,
70 * and assume that its process grid has dimension p x q.
71 * LOCr( K ) denotes the number of elements of K that a process
72 * would receive if K were distributed over the p processes of its
73 * process column.
74 * Similarly, LOCc( K ) denotes the number of elements of K that a
75 * process would receive if K were distributed over the q processes of
76 * its process row.
77 * The values of LOCr() and LOCc() may be determined via a call to the
78 * ScaLAPACK tool function, NUMROC:
79 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
80 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
81 * An upper bound for these quantities may be computed by:
82 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
83 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
84 *
85 * Arguments
86 * =========
87 *
88 * TRANS (global input) CHARACTER*1
89 * Specifies the form of the system of equations:
90 * = 'N': sub( A )*sub( X ) = sub( B )
91 * = 'T': sub( A' )*sub( X )= sub( B ), where A' is the
92 * transpose of sub( A ).
93 * = 'C': sub( A' )*sub( X )= B, where A' is the conjugate
94 * transpose of sub( A ).
95 *
96 * M (global input) INTEGER
97 * The number of rows to be operated on, i.e. the number of rows
98 * of the distributed submatrix sub( A ). M >= 0.
99 *
100 * N (global input) INTEGER
101 * The number of columns to be operated on, i.e. the number of
102 * columns of the distributed submatrix sub( A ). N >= 0.
103 *
104 * NRHS (global input) INTEGER
105 * The number of right hand sides, i.e., the number of columns
106 * of the distributed submatrix sub( B ). NRHS >= 0.
107 *
108 * A (local input) COMPLEX pointer into the local
109 * memory to an array of dimension (LLD_A,LOCc(JA+N-1)).
110 * The original M x N matrix A.
111 *
112 * IA (global input) INTEGER
113 * The row index in the global array A indicating the first
114 * row of sub( A ).
115 *
116 * JA (global input) INTEGER
117 * The column index in the global array A indicating the
118 * first column of sub( A ).
119 *
120 * DESCA (global and local input) INTEGER array of dimension DLEN_.
121 * The array descriptor for the distributed matrix A.
122 *
123 * X (local input) COMPLEX pointer into the local
124 * memory to an array of dimension (LLD_X,LOCc(JX+NRHS-1)). This
125 * array contains the local pieces of the computed solution
126 * distributed vectors for the system of linear equations.
127 *
128 * IX (global input) INTEGER
129 * The row index in the global array X indicating the first
130 * row of sub( X ).
131 *
132 * JX (global input) INTEGER
133 * The column index in the global array X indicating the
134 * first column of sub( X ).
135 *
136 * DESCX (global and local input) INTEGER array of dimension DLEN_.
137 * The array descriptor for the distributed matrix X.
138 *
139 * B (local input/local output) COMPLEX pointer into
140 * the local memory to an array of dimension
141 * (LLD_B,LOCc(JB+NRHS-1)). On entry, this array contains the
142 * local pieces of the distributes right hand side vectors for
143 * the system of linear equations. On exit, sub( B ) is over-
144 * written with the difference sub( B ) - sub( A )*sub( X ) or
145 * sub( B ) - sub( A )'*sub( X ).
146 *
147 * IB (global input) INTEGER
148 * The row index in the global array B indicating the first
149 * row of sub( B ).
150 *
151 * JB (global input) INTEGER
152 * The column index in the global array B indicating the
153 * first column of sub( B ).
154 *
155 * DESCB (global and local input) INTEGER array of dimension DLEN_.
156 * The array descriptor for the distributed matrix B.
157 *
158 * RWORK (local workspace) REAL array, dimension (LRWORK)
159 * LWORK >= Nq0 if TRANS = 'N', and LRWORK >= Mp0 otherwise.
160 *
161 * where
162 *
163 * IROFFA = MOD( IA-1, MB_A ), ICOFFA = MOD( JA-1, NB_A ),
164 * IAROW = INDXG2P( IA, MB_A, MYROW, RSRC_A, NPROW ),
165 * IACOL = INDXG2P( JA, NB_A, MYCOL, CSRC_A, NPCOL ),
166 * Mp0 = NUMROC( M+IROFFA, MB_A, MYROW, IAROW, NPROW ),
167 * Nq0 = NUMROC( N+ICOFFA, NB_A, MYCOL, IACOL, NPCOL ),
168 *
169 * INDXG2P and NUMROC are ScaLAPACK tool functions; MYROW,
170 * MYCOL, NPROW and NPCOL can be determined by calling the
171 * subroutine BLACS_GRIDINFO.
172 *
173 * RESID (global output) REAL
174 * The maximum over the number of right hand sides of
175 * norm( sub( B )- sub( A )*sub( X ) ) /
176 * ( max(m,n) * norm( sub( A ) ) * norm( sub( X ) ) * EPS ).
177 *
178 * =====================================================================
179 *
180 * .. Parameters ..
181  INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
182  \$ lld_, mb_, m_, nb_, n_, rsrc_
183  parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
184  \$ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
185  \$ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
186  REAL ZERO, ONE
187  parameter( zero = 0.0e+0, one = 1.0e+0 )
188  COMPLEX CONE
189  parameter( cone = ( 1.0e+0, 0.0e+0 ) )
190 * ..
191 * .. Local Scalars ..
192  INTEGER ICTXT, IDUMM, J, MYCOL, MYROW, N1, N2, NPCOL,
193  \$ nprow
194  REAL ANORM, BNORM, EPS, XNORM
195 * ..
196 * .. Local Arrays ..
197  REAL TEMP( 2 )
198 * ..
199 * .. External Functions ..
200  LOGICAL LSAME
201  REAL PCLANGE, PSLAMCH
202  EXTERNAL lsame, pclange, pslamch
203 * ..
204 * .. External Subroutines ..
205  EXTERNAL blacs_gridinfo, pcgemm, pscasum,
206  \$ sgamx2d
207 * ..
208 * .. Intrinsic Functions ..
209  INTRINSIC max
210 * ..
211 * .. Executable Statements ..
212 *
213 * Get grid parameters
214 *
215  ictxt = desca( ctxt_ )
216  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
217 *
218 * Quick exit if M = 0 or N = 0 or NRHS = 0
219 *
220  IF( m.LE.0 .OR. n.LE.0 .OR. nrhs.EQ.0 ) THEN
221  resid = zero
222  RETURN
223  END IF
224 *
225  IF( lsame( trans, 'T' ) .OR. lsame( trans, 'C' ) ) THEN
226  anorm = pclange( 'I', m, n, a, ia, ja, desca, rwork )
227  n1 = n
228  n2 = m
229  ELSE
230  anorm = pclange( '1', m, n, a, ia, ja, desca, rwork )
231  n1 = m
232  n2 = n
233  END IF
234 *
235  eps = pslamch( ictxt, 'Epsilon' )
236 *
237 * Compute B - sub( A )*sub( X ) (or B - sub( A' )*sub( X ) ) and
238 * store in B.
239 *
240  CALL pcgemm( trans, 'No transpose', n1, nrhs, n2, -cone, a, ia,
241  \$ ja, desca, x, ix, jx, descx, cone, b, ib, jb, descb )
242 *
243 * Compute the maximum over the number of right hand sides of
244 * norm( sub( B ) - sub( A )*sub( X ) ) /
245 * ( max(m,n) * norm( sub( A ) ) * norm( sub( X ) ) * EPS ).
246 *
247  resid = zero
248  DO 10 j = 1, nrhs
249 *
250  CALL pscasum( n1, bnorm, b, ib, jb+j-1, descb, 1 )
251  CALL pscasum( n2, xnorm, x, ix, jx+j-1, descx, 1 )
252 *
253 * Only the process columns owning the vector operands will have
254 * the correct result, the other will have zero.
255 *
256  temp( 1 ) = bnorm
257  temp( 2 ) = xnorm
258  idumm = 0
259  CALL sgamx2d( ictxt, 'All', ' ', 2, 1, temp, 2, idumm, idumm,
260  \$ -1, -1, idumm )
261  bnorm = temp( 1 )
262  xnorm = temp( 2 )
263 *
264 * Every processes have ANORM, BNORM and XNORM now.
265 *
266  IF( anorm.EQ.zero .AND. bnorm.EQ.zero ) THEN
267  resid = zero
268  ELSE IF( anorm.LE.zero .OR. xnorm.LE.zero ) THEN
269  resid = one / eps
270  ELSE
271  resid = max( resid, ( ( bnorm / anorm ) / xnorm ) /
272  \$ ( max( m, n )*eps ) )
273  END IF
274 *
275  10 CONTINUE
276 *
277  RETURN
278 *
279 * End of PCQRT16
280 *
281  END
max
#define max(A, B)
Definition: pcgemr.c:180
pcqrt16
subroutine pcqrt16(TRANS, M, N, NRHS, A, IA, JA, DESCA, X, IX, JX, DESCX, B, IB, JB, DESCB, RWORK, RESID)
Definition: pcqrt16.f:3