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__
pcagemv_(F_CHAR_T TRANS,int * M,int * N,float * ALPHA,float * A,int * IA,int * JA,int * DESCA,float * X,int * IX,int * JX,int * DESCX,int * INCX,float * BETA,float * Y,int * IY,int * JY,int * DESCY,int * INCY)20 void pcagemv_( F_CHAR_T TRANS, int * M, int * N, float * ALPHA,
21 float * A, int * IA, int * JA, int * DESCA,
22 float * X, int * IX, int * JX, int * DESCX, int * INCX,
23 float * BETA,
24 float * Y, int * IY, int * JY, int * DESCY, int * INCY )
25 #else
26 void pcagemv_( TRANS, M, N, ALPHA, A, IA, JA, DESCA, X, IX, JX, DESCX,
27 INCX, BETA, Y, IY, JY, DESCY, INCY )
28 /*
29 * .. Scalar Arguments ..
30 */
31 F_CHAR_T TRANS;
32 int * IA, * INCX, * INCY, * IX, * IY, * JA, * JX, * JY,
33 * M, * N;
34 float * ALPHA, * BETA;
35 /*
36 * .. Array Arguments ..
37 */
38 int * DESCA, * DESCX, * DESCY;
39 float * A, * X, * Y;
40 #endif
41 {
42 /*
43 * Purpose
44 * =======
45 *
46 * PCAGEMV performs one of the matrix-vector operations
47 *
48 * sub( Y ) := abs( alpha )*abs( sub( A ) )*abs( sub( X ) ) +
49 * abs( beta*sub( Y ) ),
50 * or
51 *
52 * sub( Y ) := abs( alpha )*abs( sub( A )' )*abs( sub( X ) ) +
53 * abs( beta*sub( Y ) ),
54 * or
55 *
56 * sub( Y ) := abs( alpha )*abs( conjg( sub( A )' ) )*abs( sub( X ) )
57 * + abs( beta*sub( Y ) ),
58 *
59 * where
60 *
61 * sub( A ) denotes A(IA:IA+M-1,JA:JA+N-1).
62 *
63 * When TRANS = 'N',
64 *
65 * sub( X ) denotes X(IX:IX,JX:JX+N-1), if INCX = M_X,
66 * X(IX:IX+N-1,JX:JX), if INCX = 1 and INCX <> M_X,
67 * and,
68 *
69 * sub( Y ) denotes Y(IY:IY,JY:JY+M-1), if INCY = M_Y,
70 * Y(IY:IY+M-1,JY:JY), if INCY = 1 and INCY <> M_Y,
71 * and, otherwise
72 *
73 * sub( X ) denotes X(IX:IX,JX:JX+M-1), if INCX = M_X,
74 * X(IX:IX+M-1,JX:JX), if INCX = 1 and INCX <> M_X,
75 * and,
76 *
77 * sub( Y ) denotes Y(IY:IY,JY:JY+N-1), if INCY = M_Y,
78 * Y(IY:IY+N-1,JY:JY), if INCY = 1 and INCY <> M_Y.
79 *
80 * Alpha and beta are real scalars, sub( Y ) is a real subvector,
81 * sub( X ) is a subvector and sub( A ) is an m by n submatrix.
82 *
83 * Notes
84 * =====
85 *
86 * A description vector is associated with each 2D block-cyclicly dis-
87 * tributed matrix. This vector stores the information required to
88 * establish the mapping between a matrix entry and its corresponding
89 * process and memory location.
90 *
91 * In the following comments, the character _ should be read as
92 * "of the distributed matrix". Let A be a generic term for any 2D
93 * block cyclicly distributed matrix. Its description vector is DESC_A:
94 *
95 * NOTATION STORED IN EXPLANATION
96 * ---------------- --------------- ------------------------------------
97 * DTYPE_A (global) DESCA[ DTYPE_ ] The descriptor type.
98 * CTXT_A (global) DESCA[ CTXT_ ] The BLACS context handle, indicating
99 * the NPROW x NPCOL BLACS process grid
100 * A is distributed over. The context
101 * itself is global, but the handle
102 * (the integer value) may vary.
103 * M_A (global) DESCA[ M_ ] The number of rows in the distribu-
104 * ted matrix A, M_A >= 0.
105 * N_A (global) DESCA[ N_ ] The number of columns in the distri-
106 * buted matrix A, N_A >= 0.
107 * IMB_A (global) DESCA[ IMB_ ] The number of rows of the upper left
108 * block of the matrix A, IMB_A > 0.
109 * INB_A (global) DESCA[ INB_ ] The number of columns of the upper
110 * left block of the matrix A,
111 * INB_A > 0.
112 * MB_A (global) DESCA[ MB_ ] The blocking factor used to distri-
113 * bute the last M_A-IMB_A rows of A,
114 * MB_A > 0.
115 * NB_A (global) DESCA[ NB_ ] The blocking factor used to distri-
116 * bute the last N_A-INB_A columns of
117 * A, NB_A > 0.
118 * RSRC_A (global) DESCA[ RSRC_ ] The process row over which the first
119 * row of the matrix A is distributed,
120 * NPROW > RSRC_A >= 0.
121 * CSRC_A (global) DESCA[ CSRC_ ] The process column over which the
122 * first column of A is distributed.
123 * NPCOL > CSRC_A >= 0.
124 * LLD_A (local) DESCA[ LLD_ ] The leading dimension of the local
125 * array storing the local blocks of
126 * the distributed matrix A,
127 * IF( Lc( 1, N_A ) > 0 )
128 * LLD_A >= MAX( 1, Lr( 1, M_A ) )
129 * ELSE
130 * LLD_A >= 1.
131 *
132 * Let K be the number of rows of a matrix A starting at the global in-
133 * dex IA,i.e, A( IA:IA+K-1, : ). Lr( IA, K ) denotes the number of rows
134 * that the process of row coordinate MYROW ( 0 <= MYROW < NPROW ) would
135 * receive if these K rows were distributed over NPROW processes. If K
136 * is the number of columns of a matrix A starting at the global index
137 * JA, i.e, A( :, JA:JA+K-1, : ), Lc( JA, K ) denotes the number of co-
138 * lumns that the process MYCOL ( 0 <= MYCOL < NPCOL ) would receive if
139 * these K columns were distributed over NPCOL processes.
140 *
141 * The values of Lr() and Lc() may be determined via a call to the func-
142 * tion PB_Cnumroc:
143 * Lr( IA, K ) = PB_Cnumroc( K, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW )
144 * Lc( JA, K ) = PB_Cnumroc( K, JA, INB_A, NB_A, MYCOL, CSRC_A, NPCOL )
145 *
146 * Arguments
147 * =========
148 *
149 * TRANS (global input) CHARACTER*1
150 * On entry, TRANS specifies the operation to be performed as
151 * follows:
152 *
153 * TRANS = 'N' or 'n'
154 * sub( Y ) := |alpha|*|sub( A ) | * |sub( X )| +
155 * |beta*sub( Y )|,
156 *
157 * TRANS = 'T' or 't',
158 * sub( Y ) := |alpha|*|sub( A )'| * |sub( X )| +
159 * |beta*sub( Y )|,
160 *
161 * TRANS = 'C' or 'c',
162 * sub( Y ) := |alpha|*|conjg( sub( A )' )|*|sub( X )| +
163 * |beta*sub( Y )|.
164 *
165 * M (global input) INTEGER
166 * On entry, M specifies the number of rows of the submatrix
167 * sub( A ). M must be at least zero.
168 *
169 * N (global input) INTEGER
170 * On entry, N specifies the number of columns of the submatrix
171 * sub( A ). N must be at least zero.
172 *
173 * ALPHA (global input) REAL
174 * On entry, ALPHA specifies the scalar alpha. When ALPHA is
175 * supplied as zero then the local entries of the arrays A
176 * and X corresponding to the entries of the submatrix sub( A )
177 * and the subvector sub( X ) need not be set on input.
178 *
179 * A (local input) COMPLEX array
180 * On entry, A is an array of dimension (LLD_A, Ka), where Ka is
181 * at least Lc( 1, JA+N-1 ). Before entry, this array contains
182 * the local entries of the matrix A.
183 *
184 * IA (global input) INTEGER
185 * On entry, IA specifies A's global row index, which points to
186 * the beginning of the submatrix sub( A ).
187 *
188 * JA (global input) INTEGER
189 * On entry, JA specifies A's global column index, which points
190 * to the beginning of the submatrix sub( A ).
191 *
192 * DESCA (global and local input) INTEGER array
193 * On entry, DESCA is an integer array of dimension DLEN_. This
194 * is the array descriptor for the matrix A.
195 *
196 * X (local input) COMPLEX array
197 * On entry, X is an array of dimension (LLD_X, Kx), where LLD_X
198 * is at least MAX( 1, Lr( 1, IX ) ) when INCX = M_X and
199 * MAX( 1, Lr( 1, IX+Lx-1 ) ) otherwise, and, Kx is at least
200 * Lc( 1, JX+Lx-1 ) when INCX = M_X and Lc( 1, JX ) otherwise.
201 * Lx is N when TRANS = 'N' or 'n' and M otherwise. Before en-
202 * try, this array contains the local entries of the matrix X.
203 *
204 * IX (global input) INTEGER
205 * On entry, IX specifies X's global row index, which points to
206 * the beginning of the submatrix sub( X ).
207 *
208 * JX (global input) INTEGER
209 * On entry, JX specifies X's global column index, which points
210 * to the beginning of the submatrix sub( X ).
211 *
212 * DESCX (global and local input) INTEGER array
213 * On entry, DESCX is an integer array of dimension DLEN_. This
214 * is the array descriptor for the matrix X.
215 *
216 * INCX (global input) INTEGER
217 * On entry, INCX specifies the global increment for the
218 * elements of X. Only two values of INCX are supported in
219 * this version, namely 1 and M_X. INCX must not be zero.
220 *
221 * BETA (global input) REAL
222 * On entry, BETA specifies the scalar beta. When BETA is
223 * supplied as zero then the local entries of the array Y
224 * corresponding to the entries of the subvector sub( Y ) need
225 * not be set on input.
226 *
227 * Y (local input/local output) REAL array
228 * On entry, Y is an array of dimension (LLD_Y, Ky), where LLD_Y
229 * is at least MAX( 1, Lr( 1, IY ) ) when INCY = M_Y and
230 * MAX( 1, Lr( 1, IY+Ly-1 ) ) otherwise, and, Ky is at least
231 * Lc( 1, JY+Ly-1 ) when INCY = M_Y and Lc( 1, JY ) otherwise.
232 * Ly is M when TRANS = 'N' or 'n' and N otherwise. Before en-
233 * try, this array contains the local entries of the matrix Y.
234 * On exit, sub( Y ) is overwritten by the updated subvector.
235 *
236 * IY (global input) INTEGER
237 * On entry, IY specifies Y's global row index, which points to
238 * the beginning of the submatrix sub( Y ).
239 *
240 * JY (global input) INTEGER
241 * On entry, JY specifies Y's global column index, which points
242 * to the beginning of the submatrix sub( Y ).
243 *
244 * DESCY (global and local input) INTEGER array
245 * On entry, DESCY is an integer array of dimension DLEN_. This
246 * is the array descriptor for the matrix Y.
247 *
248 * INCY (global input) INTEGER
249 * On entry, INCY specifies the global increment for the
250 * elements of Y. Only two values of INCY are supported in
251 * this version, namely 1 and M_Y. INCY must not be zero.
252 *
253 * -- Written on April 1, 1998 by
254 * Antoine Petitet, University of Tennessee, Knoxville 37996, USA.
255 *
256 * ---------------------------------------------------------------------
257 */
258 /*
259 * .. Local Scalars ..
260 */
261 char TrA, Yroc, * one, * tbeta, top;
262 int Acol, Ai, Aii, Aimb1, Ainb1, Aj, Ajj, Ald, Amb, Amp, Anb,
263 Anq, Arow, XAfr, Xi, Xj, YAfr, YApbY, YAsum, Ycol, Yi, Yii,
264 Yj, Yjj, Yld, Ynp, Ynq, Yrow, ctxt, info, ione=1, mycol,
265 myrow, nota, npcol, nprow;
266 PBTYP_T * type, * utyp;
267 /*
268 * .. Local Arrays ..
269 */
270 int Ad [DLEN_], Ad0[DLEN_], XAd[DLEN_], Xd[DLEN_], YAd[DLEN_],
271 Yd [DLEN_];
272 char * XA = NULL, * YA = NULL;
273 /* ..
274 * .. Executable Statements ..
275 *
276 */
277 nota = ( ( TrA = Mupcase( F2C_CHAR( TRANS )[0] ) ) == CNOTRAN );
278 PB_CargFtoC( *IA, *JA, DESCA, &Ai, &Aj, Ad );
279 PB_CargFtoC( *IX, *JX, DESCX, &Xi, &Xj, Xd );
280 PB_CargFtoC( *IY, *JY, DESCY, &Yi, &Yj, Yd );
281 #ifndef NO_ARGCHK
282 /*
283 * Test the input parameters
284 */
285 Cblacs_gridinfo( ( ctxt = Ad[CTXT_] ), &nprow, &npcol, &myrow, &mycol );
286 if( !( info = ( ( nprow == -1 ) ? -( 801 + CTXT_ ) : 0 ) ) )
287 {
288 if( ( !nota ) && ( TrA != CTRAN ) && ( TrA != CCOTRAN ) )
289 {
290 PB_Cwarn( ctxt, __LINE__, "PCAGEMV", "Illegal TRANS=%c\n", TrA );
291 info = -1;
292 }
293 PB_Cchkmat( ctxt, "PCAGEMV", "A", *M, 2, *N, 3, Ai, Aj, Ad, 8,
294 &info );
295 if( nota )
296 {
297 PB_Cchkvec( ctxt, "PCAGEMV", "X", *N, 3, Xi, Xj, Xd, *INCX, 12,
298 &info );
299 PB_Cchkvec( ctxt, "PCAGEMV", "Y", *M, 2, Yi, Yj, Yd, *INCY, 18,
300 &info );
301 }
302 else
303 {
304 PB_Cchkvec( ctxt, "PCAGEMV", "X", *M, 2, Xi, Xj, Xd, *INCX, 12,
305 &info );
306 PB_Cchkvec( ctxt, "PCAGEMV", "Y", *N, 3, Yi, Yj, Yd, *INCY, 18,
307 &info );
308 }
309 }
310 if( info ) { PB_Cabort( ctxt, "PCAGEMV", info ); return; }
311 #endif
312 /*
313 * Quick return if possible
314 */
315 if( ( *M == 0 ) || ( *N == 0 ) ||
316 ( ( ALPHA[REAL_PART] == ZERO ) && ( BETA[REAL_PART] == ONE ) ) )
317 return;
318 /*
319 * Retrieve process grid information
320 */
321 #ifdef NO_ARGCHK
322 Cblacs_gridinfo( ( ctxt = Ad[CTXT_] ), &nprow, &npcol, &myrow, &mycol );
323 #endif
324 /*
325 * Get type structure
326 */
327 type = PB_Cctypeset();
328 utyp = PB_Cstypeset();
329 /*
330 * When alpha is zero
331 */
332 if( ALPHA[REAL_PART] == ZERO )
333 {
334 /*
335 * Retrieve sub( Y )'s local information: Yii, Yjj, Yrow, Ycol
336 */
337 PB_Cinfog2l( Yi, Yj, Yd, nprow, npcol, myrow, mycol, &Yii, &Yjj,
338 &Yrow, &Ycol );
339
340 if( *INCY == Yd[M_] )
341 {
342 /*
343 * sub( Y ) resides in (a) process row(s)
344 */
345 if( ( myrow == Yrow ) || ( Yrow < 0 ) )
346 {
347 /*
348 * Make sure I own some data and scale sub( Y )
349 */
350 Ynq = PB_Cnumroc( ( nota ? *M : *N ), Yj, Yd[INB_], Yd[NB_], mycol,
351 Yd[CSRC_], npcol );
352 if( Ynq > 0 )
353 {
354 Yld = Yd[LLD_];
355 sascal_( &Ynq, ((char *) BETA), Mptr( ((char *) Y), Yii,
356 Yjj, Yld, utyp->size ), &Yld );
357 }
358 }
359 }
360 else
361 {
362 /*
363 * sub( Y ) resides in (a) process column(s)
364 */
365 if( ( mycol == Ycol ) || ( Ycol < 0 ) )
366 {
367 /*
368 * Make sure I own some data and scale sub( Y )
369 */
370 Ynp = PB_Cnumroc( ( nota ? *M : *N ), Yi, Yd[IMB_], Yd[MB_], myrow,
371 Yd[RSRC_], nprow );
372 if( Ynp > 0 )
373 {
374 sascal_( &Ynp, ((char *) BETA), Mptr( ((char *) Y), Yii,
375 Yjj, Yd[LLD_], utyp->size ), INCY );
376 }
377 }
378 }
379 return;
380 }
381 /*
382 * Compute descriptor Ad0 for sub( A )
383 */
384 PB_Cdescribe( *M, *N, Ai, Aj, Ad, nprow, npcol, myrow, mycol, &Aii, &Ajj,
385 &Ald, &Aimb1, &Ainb1, &Amb, &Anb, &Arow, &Acol, Ad0 );
386
387 Yroc = ( *INCY == Yd[M_] ? CROW : CCOLUMN );
388
389 if( nota )
390 {
391 /*
392 * Reuse sub( Y ) and/or create vector YA in process columns spanned by sub( A )
393 */
394 PB_CInOutV( utyp, COLUMN, *M, *N, Ad0, 1, ((char *) BETA), ((char *) Y),
395 Yi, Yj, Yd, &Yroc, &tbeta, &YA, YAd, &YAfr, &YAsum, &YApbY );
396 /*
397 * Replicate sub( X ) in process rows spanned by sub( A ) -> XA
398 */
399 PB_CInV( type, NOCONJG, ROW, *M, *N, Ad0, 1, ((char *) X), Xi, Xj, Xd,
400 ( *INCX == Xd[M_] ? ROW : COLUMN ), &XA, XAd, &XAfr );
401 /*
402 * Local matrix-vector multiply iff I own some data
403 */
404 Amp = PB_Cnumroc( *M, 0, Ad0[IMB_], Ad0[MB_], myrow, Ad0[RSRC_], nprow );
405 Anq = PB_Cnumroc( *N, 0, Ad0[INB_], Ad0[NB_], mycol, Ad0[CSRC_], npcol );
406 if( ( Amp > 0 ) && ( Anq > 0 ) )
407 {
408 cagemv_( TRANS, &Amp, &Anq, ((char *) ALPHA), Mptr( ((char *) A),
409 Aii, Ajj, Ald, type->size), &Ald, XA, &XAd[LLD_], tbeta,
410 YA, &ione );
411 }
412 if( XAfr ) free( XA );
413 /*
414 * Combine the partial column results into YA
415 */
416 if( YAsum && ( Amp > 0 ) )
417 {
418 top = *PB_Ctop( &ctxt, COMBINE, ROW, TOP_GET );
419 Csgsum2d( ctxt, ROW, &top, Amp, 1, YA, YAd[LLD_], myrow,
420 YAd[CSRC_] );
421 }
422 }
423 else
424 {
425 /*
426 * Reuse sub( Y ) and/or create vector YA in process rows spanned by sub( A )
427 */
428 PB_CInOutV( utyp, ROW, *M, *N, Ad0, 1, ((char *) BETA), ((char *) Y), Yi,
429 Yj, Yd, &Yroc, &tbeta, &YA, YAd, &YAfr, &YAsum, &YApbY );
430 /*
431 * Replicate sub( X ) in process columns spanned by sub( A ) -> XA
432 */
433 PB_CInV( type, NOCONJG, COLUMN, *M, *N, Ad0, 1, ((char *) X), Xi, Xj, Xd,
434 ( *INCX == Xd[M_] ? ROW : COLUMN ), &XA, XAd, &XAfr );
435 /*
436 * Local matrix-vector multiply iff I own some data
437 */
438 Amp = PB_Cnumroc( *M, 0, Ad0[IMB_], Ad0[MB_], myrow, Ad0[RSRC_], nprow );
439 Anq = PB_Cnumroc( *N, 0, Ad0[INB_], Ad0[NB_], mycol, Ad0[CSRC_], npcol );
440 if( ( Amp > 0 ) && ( Anq > 0 ) )
441 {
442 cagemv_( TRANS, &Amp, &Anq, ((char *) ALPHA), Mptr( ((char *) A),
443 Aii, Ajj, Ald, type->size ), &Ald, XA, &ione, tbeta, YA,
444 &YAd[LLD_] );
445 }
446 if( XAfr ) free( XA );
447 /*
448 * Combine the partial row results into YA
449 */
450 if( YAsum && ( Anq > 0 ) )
451 {
452 top = *PB_Ctop( &ctxt, COMBINE, COLUMN, TOP_GET );
453 Csgsum2d( ctxt, COLUMN, &top, 1, Anq, YA, YAd[LLD_], YAd[RSRC_],
454 mycol );
455 }
456 }
457 /*
458 * sub( Y ) := beta * sub( Y ) + YA (if necessary)
459 */
460 if( YApbY )
461 {
462 /*
463 * Retrieve sub( Y )'s local information: Yii, Yjj, Yrow, Ycol
464 */
465 PB_Cinfog2l( Yi, Yj, Yd, nprow, npcol, myrow, mycol, &Yii, &Yjj, &Yrow,
466 &Ycol );
467
468 if( *INCY == Yd[M_] )
469 {
470 /*
471 * sub( Y ) resides in (a) process row(s)
472 */
473 if( ( myrow == Yrow ) || ( Yrow < 0 ) )
474 {
475 /*
476 * Make sure I own some data and scale sub( Y )
477 */
478 Ynq = PB_Cnumroc( ( nota ? *M : *N ), Yj, Yd[INB_], Yd[NB_], mycol,
479 Yd[CSRC_], npcol );
480 if( Ynq > 0 )
481 {
482 Yld = Yd[LLD_];
483 sascal_( &Ynq, ((char *) BETA), Mptr( ((char *) Y), Yii,
484 Yjj, Yld, utyp->size ), &Yld );
485 }
486 }
487 }
488 else
489 {
490 /*
491 * sub( Y ) resides in (a) process column(s)
492 */
493 if( ( mycol == Ycol ) || ( Ycol < 0 ) )
494 {
495 /*
496 * Make sure I own some data and scale sub( Y )
497 */
498 Ynp = PB_Cnumroc( ( nota ? *M : *N ), Yi, Yd[IMB_], Yd[MB_], myrow,
499 Yd[RSRC_], nprow );
500 if( Ynp > 0 )
501 {
502 sascal_( &Ynp, ((char *) BETA), Mptr( ((char *) Y), Yii,
503 Yjj, Yd[LLD_], utyp->size ), INCY );
504 }
505 }
506 }
507
508 one = utyp->one;
509
510 if( nota )
511 {
512 PB_Cpaxpby( utyp, NOCONJG, *M, 1, one, YA, 0, 0, YAd, COLUMN, one,
513 ((char *) Y), Yi, Yj, Yd, &Yroc );
514 }
515 else
516 {
517 PB_Cpaxpby( utyp, NOCONJG, 1, *N, one, YA, 0, 0, YAd, ROW, one,
518 ((char *) Y), Yi, Yj, Yd, &Yroc );
519 }
520 }
521 if( YAfr ) free( YA );
522 /*
523 * End of PCAGEMV
524 */
525 }
526