1 /* ---------------------------------------------------------------------
2 *
3 * -- PBLAS routine (version 2.0) --
4 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
5 * and University of California, Berkeley.
6 * April 1, 1998
7 *
8 * ---------------------------------------------------------------------
9 */
10 /*
11 * Include files
12 */
13 #include "pblas.h"
14 #include "PBpblas.h"
15 #include "PBtools.h"
16 #include "PBblacs.h"
17 #include "PBblas.h"
18
19 #ifdef __STDC__
psasum_(int * N,float * ASUM,float * X,int * IX,int * JX,int * DESCX,int * INCX)20 void psasum_( int * N, float * ASUM,
21 float * X, int * IX, int * JX, int * DESCX, int * INCX )
22 #else
23 void psasum_( N, ASUM, X, IX, JX, DESCX, INCX )
24 /*
25 * .. Scalar Arguments ..
26 */
27 int * INCX, * IX, * JX, * N;
28 float * ASUM;
29 /*
30 * .. Array Arguments ..
31 */
32 int * DESCX;
33 float * X;
34 #endif
35 {
36 /*
37 * Purpose
38 * =======
39 *
40 * PSASUM returns the sum of absolute values of the entries of a subvec-
41 * tor sub( X ),
42 *
43 * where
44 *
45 * sub( X ) denotes X(IX,JX:JX+N-1) if INCX = M_X,
46 * X(IX:IX+N-1,JX) if INCX = 1 and INCX <> M_X.
47 *
48 * Notes
49 * =====
50 *
51 * A description vector is associated with each 2D block-cyclicly dis-
52 * tributed matrix. This vector stores the information required to
53 * establish the mapping between a matrix entry and its corresponding
54 * process and memory location.
55 *
56 * In the following comments, the character _ should be read as
57 * "of the distributed matrix". Let A be a generic term for any 2D
58 * block cyclicly distributed matrix. Its description vector is DESC_A:
59 *
60 * NOTATION STORED IN EXPLANATION
61 * ---------------- --------------- ------------------------------------
62 * DTYPE_A (global) DESCA[ DTYPE_ ] The descriptor type.
63 * CTXT_A (global) DESCA[ CTXT_ ] The BLACS context handle, indicating
64 * the NPROW x NPCOL BLACS process grid
65 * A is distributed over. The context
66 * itself is global, but the handle
67 * (the integer value) may vary.
68 * M_A (global) DESCA[ M_ ] The number of rows in the distribu-
69 * ted matrix A, M_A >= 0.
70 * N_A (global) DESCA[ N_ ] The number of columns in the distri-
71 * buted matrix A, N_A >= 0.
72 * IMB_A (global) DESCA[ IMB_ ] The number of rows of the upper left
73 * block of the matrix A, IMB_A > 0.
74 * INB_A (global) DESCA[ INB_ ] The number of columns of the upper
75 * left block of the matrix A,
76 * INB_A > 0.
77 * MB_A (global) DESCA[ MB_ ] The blocking factor used to distri-
78 * bute the last M_A-IMB_A rows of A,
79 * MB_A > 0.
80 * NB_A (global) DESCA[ NB_ ] The blocking factor used to distri-
81 * bute the last N_A-INB_A columns of
82 * A, NB_A > 0.
83 * RSRC_A (global) DESCA[ RSRC_ ] The process row over which the first
84 * row of the matrix A is distributed,
85 * NPROW > RSRC_A >= 0.
86 * CSRC_A (global) DESCA[ CSRC_ ] The process column over which the
87 * first column of A is distributed.
88 * NPCOL > CSRC_A >= 0.
89 * LLD_A (local) DESCA[ LLD_ ] The leading dimension of the local
90 * array storing the local blocks of
91 * the distributed matrix A,
92 * IF( Lc( 1, N_A ) > 0 )
93 * LLD_A >= MAX( 1, Lr( 1, M_A ) )
94 * ELSE
95 * LLD_A >= 1.
96 *
97 * Let K be the number of rows of a matrix A starting at the global in-
98 * dex IA,i.e, A( IA:IA+K-1, : ). Lr( IA, K ) denotes the number of rows
99 * that the process of row coordinate MYROW ( 0 <= MYROW < NPROW ) would
100 * receive if these K rows were distributed over NPROW processes. If K
101 * is the number of columns of a matrix A starting at the global index
102 * JA, i.e, A( :, JA:JA+K-1, : ), Lc( JA, K ) denotes the number of co-
103 * lumns that the process MYCOL ( 0 <= MYCOL < NPCOL ) would receive if
104 * these K columns were distributed over NPCOL processes.
105 *
106 * The values of Lr() and Lc() may be determined via a call to the func-
107 * tion PB_Cnumroc:
108 * Lr( IA, K ) = PB_Cnumroc( K, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW )
109 * Lc( JA, K ) = PB_Cnumroc( K, JA, INB_A, NB_A, MYCOL, CSRC_A, NPCOL )
110 *
111 * Arguments
112 * =========
113 *
114 * N (global input) INTEGER
115 * On entry, N specifies the length of the subvector sub( X ).
116 * N must be at least zero.
117 *
118 * ASUM (local output) REAL
119 * On exit, ASUM specifies the sum of absolute values of the
120 * subvector sub( X ) only in its scope (See below for further
121 * details).
122 *
123 * X (local input) REAL array
124 * On entry, X is an array of dimension (LLD_X, Kx), where LLD_X
125 * is at least MAX( 1, Lr( 1, IX ) ) when INCX = M_X and
126 * MAX( 1, Lr( 1, IX+N-1 ) ) otherwise, and, Kx is at least
127 * Lc( 1, JX+N-1 ) when INCX = M_X and Lc( 1, JX ) otherwise.
128 * Before entry, this array contains the local entries of the
129 * matrix X.
130 *
131 * IX (global input) INTEGER
132 * On entry, IX specifies X's global row index, which points to
133 * the beginning of the submatrix sub( X ).
134 *
135 * JX (global input) INTEGER
136 * On entry, JX specifies X's global column index, which points
137 * to the beginning of the submatrix sub( X ).
138 *
139 * DESCX (global and local input) INTEGER array
140 * On entry, DESCX is an integer array of dimension DLEN_. This
141 * is the array descriptor for the matrix X.
142 *
143 * INCX (global input) INTEGER
144 * On entry, INCX specifies the global increment for the
145 * elements of X. Only two values of INCX are supported in
146 * this version, namely 1 and M_X. INCX must not be zero.
147 *
148 * Further Details
149 * ===============
150 *
151 * When the result of a vector-oriented PBLAS call is a scalar, this
152 * scalar is set only within the process scope which owns the vector(s)
153 * being operated on. Let sub( X ) be a generic term for the input vec-
154 * tor(s). Then, the processes owning the correct the answer is determi-
155 * ned as follows: if an operation involves more than one vector, the
156 * processes receiving the result will be the union of the following set
157 * of processes for each vector:
158 *
159 * If N = 1, M_X = 1 and INCX = 1, then one cannot determine if a pro-
160 * cess row or process column owns the vector operand, therefore only
161 * the process owning sub( X ) receives the correct result;
162 *
163 * If INCX = M_X, then sub( X ) is a vector distributed over a process
164 * row. Each process in this row receives the result;
165 *
166 * If INCX = 1, then sub( X ) is a vector distributed over a process
167 * column. Each process in this column receives the result;
168 *
169 * -- Written on April 1, 1998 by
170 * Antoine Petitet, University of Tennessee, Knoxville 37996, USA.
171 *
172 * ---------------------------------------------------------------------
173 */
174 /*
175 * .. Local Scalars ..
176 */
177 char top;
178 int Xcol, Xi, Xii, Xj, Xjj, Xld, Xnp, Xnq, Xrow, ctxt, info,
179 mycol, myrow, npcol, nprow;
180 /*
181 * .. Local Arrays ..
182 */
183 int Xd[DLEN_];
184 /* ..
185 * .. Executable Statements ..
186 *
187 */
188 PB_CargFtoC( *IX, *JX, DESCX, &Xi, &Xj, Xd );
189 #ifndef NO_ARGCHK
190 /*
191 * Test the input parameters
192 */
193 Cblacs_gridinfo( ( ctxt = Xd[CTXT_] ), &nprow, &npcol, &myrow, &mycol );
194 if( !( info = ( ( nprow == -1 ) ? -( 601 + CTXT_ ) : 0 ) ) )
195 PB_Cchkvec( ctxt, "PSASUM", "X", *N, 1, Xi, Xj, Xd, *INCX, 6, &info );
196 if( info ) { PB_Cabort( ctxt, "PSASUM", info ); return; }
197 #endif
198 /*
199 * Initialize ASUM
200 */
201 *ASUM = ZERO;
202 /*
203 * Quick return if possible
204 */
205 if( *N == 0 ) return;
206 /*
207 * Retrieve process grid information
208 */
209 #ifdef NO_ARGCHK
210 Cblacs_gridinfo( ( ctxt = Xd[CTXT_] ), &nprow, &npcol, &myrow, &mycol );
211 #endif
212 /*
213 * Retrieve sub( X )'s local information: Xii, Xjj, Xrow, Xcol
214 */
215 PB_Cinfog2l( Xi, Xj, Xd, nprow, npcol, myrow, mycol, &Xii, &Xjj,
216 &Xrow, &Xcol );
217 /*
218 * Handle degenerate case separately, sub( X )'s scope is just one process
219 */
220 if( ( *N == 1 ) && ( *INCX == 1 ) && ( Xd[M_] == 1 ) )
221 {
222 /*
223 * Make sure I own some data and compute ASUM
224 */
225 if( ( ( myrow == Xrow ) || ( Xrow < 0 ) ) &&
226 ( ( mycol == Xcol ) || ( Xcol < 0 ) ) )
227 {
228 *ASUM = ABS( X[Xii+Xjj*Xd[LLD_]] );
229 }
230 return;
231 }
232 else if( *INCX == Xd[M_] )
233 {
234 /*
235 * sub( X ) resides in (a) process row(s)
236 */
237 if( ( myrow == Xrow ) || ( Xrow < 0 ) )
238 {
239 /*
240 * Make sure I own some data and compute the local sum
241 */
242 Xnq = PB_Cnumroc( *N, Xj, Xd[INB_], Xd[NB_], mycol, Xd[CSRC_], npcol );
243 if( Xnq > 0 )
244 {
245 Xld = Xd[LLD_];
246 svasum_( &Xnq, ((char *) ASUM), ((char *)( X+(Xii+Xjj*Xld) )),
247 &Xld );
248 }
249 /*
250 * If Xnq <= 0, ASUM is zero (see initialization above)
251 */
252 if( ( npcol > 1 ) && ( Xcol >= 0 ) )
253 {
254 /*
255 * Combine the local results if npcol > 1 and Xcol >= 0, i.e sub( X ) is
256 * distributed.
257 */
258 top = *PB_Ctop( &ctxt, COMBINE, ROW, TOP_GET );
259 Csgsum2d( ctxt, ROW, &top, 1, 1, ((char *)ASUM), 1, -1,
260 mycol );
261 }
262 }
263 return;
264 }
265 else
266 {
267 /*
268 * sub( X ) resides in (a) process column(s)
269 */
270 if( ( mycol == Xcol ) || ( Xcol < 0 ) )
271 {
272 /*
273 * Make sure I own some data and compute the local sum
274 */
275 Xnp = PB_Cnumroc( *N, Xi, Xd[IMB_], Xd[MB_], myrow, Xd[RSRC_], nprow );
276 if( Xnp > 0 )
277 {
278 svasum_( &Xnp, ((char *) ASUM),
279 ((char *)( X+(Xii+Xjj*Xd[LLD_]) )), INCX );
280 }
281 /*
282 * If Xnp <= 0, ASUM is zero (see initialization above)
283 */
284 if( ( nprow > 1 ) && ( Xrow >= 0 ) )
285 {
286 /*
287 * Combine the local results if nprow > 1 and Xrow >= 0, i.e sub( X ) is
288 * distributed.
289 */
290 top = *PB_Ctop( &ctxt, COMBINE, COLUMN, TOP_GET );
291 Csgsum2d( ctxt, COLUMN, &top, 1, 1, ((char *)ASUM), 1, -1,
292 mycol );
293 }
294 }
295 return;
296 }
297 /*
298 * End of PSASUM
299 */
300 }
301