1 /* ---------------------------------------------------------------------
2 *
3 * -- PBLAS auxiliary 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__
PB_CInV(PBTYP_T * TYPE,char * CONJUG,char * ROWCOL,int M,int N,int * DESCA,int K,char * X,int IX,int JX,int * DESCX,char * XROC,char ** XAPTR,int * DXA,int * XAFREE)20 void PB_CInV( PBTYP_T * TYPE, char * CONJUG, char * ROWCOL, int M,
21 int N, int * DESCA, int K, char * X, int IX, int JX,
22 int * DESCX, char * XROC, char * * XAPTR, int * DXA,
23 int * XAFREE )
24 #else
25 void PB_CInV( TYPE, CONJUG, ROWCOL, M, N, DESCA, K, X, IX, JX, DESCX,
26 XROC, XAPTR, DXA, XAFREE )
27 /*
28 * .. Scalar Arguments ..
29 */
30 char * CONJUG, * ROWCOL, * XROC;
31 int * XAFREE, IX, JX, K, M, N;
32 PBTYP_T * TYPE;
33 /*
34 * .. Array Arguments ..
35 */
36 int * DESCA, * DESCX, * DXA;
37 char * X, * * XAPTR;
38 #endif
39 {
40 /*
41 * Purpose
42 * =======
43 *
44 * PB_CInV returns a pointer to an array that contains a one-dimensional
45 * input only subvector which is replicated over the rows or columns of
46 * a submatrix described by DESCA. A subvector is specified on input to
47 * this routine that is reused whenever possible. On return, the subvec-
48 * tor is specified by a pointer to some data, a descriptor array des-
49 * cribing its layout and a logical value indicating if this local piece
50 * of data has been dynamically allocated by this function. This routine
51 * is specifically designed for traditional Level 2 like PBLAS opera-
52 * tions using an input only vector such as PxGER, PxSYR ...
53 *
54 * Notes
55 * =====
56 *
57 * A description vector is associated with each 2D block-cyclicly dis-
58 * tributed matrix. This vector stores the information required to
59 * establish the mapping between a matrix entry and its corresponding
60 * process and memory location.
61 *
62 * In the following comments, the character _ should be read as
63 * "of the distributed matrix". Let A be a generic term for any 2D
64 * block cyclicly distributed matrix. Its description vector is DESC_A:
65 *
66 * NOTATION STORED IN EXPLANATION
67 * ---------------- --------------- ------------------------------------
68 * DTYPE_A (global) DESCA[ DTYPE_ ] The descriptor type.
69 * CTXT_A (global) DESCA[ CTXT_ ] The BLACS context handle, indicating
70 * the NPROW x NPCOL BLACS process grid
71 * A is distributed over. The context
72 * itself is global, but the handle
73 * (the integer value) may vary.
74 * M_A (global) DESCA[ M_ ] The number of rows in the distribu-
75 * ted matrix A, M_A >= 0.
76 * N_A (global) DESCA[ N_ ] The number of columns in the distri-
77 * buted matrix A, N_A >= 0.
78 * IMB_A (global) DESCA[ IMB_ ] The number of rows of the upper left
79 * block of the matrix A, IMB_A > 0.
80 * INB_A (global) DESCA[ INB_ ] The number of columns of the upper
81 * left block of the matrix A,
82 * INB_A > 0.
83 * MB_A (global) DESCA[ MB_ ] The blocking factor used to distri-
84 * bute the last M_A-IMB_A rows of A,
85 * MB_A > 0.
86 * NB_A (global) DESCA[ NB_ ] The blocking factor used to distri-
87 * bute the last N_A-INB_A columns of
88 * A, NB_A > 0.
89 * RSRC_A (global) DESCA[ RSRC_ ] The process row over which the first
90 * row of the matrix A is distributed,
91 * NPROW > RSRC_A >= 0.
92 * CSRC_A (global) DESCA[ CSRC_ ] The process column over which the
93 * first column of A is distributed.
94 * NPCOL > CSRC_A >= 0.
95 * LLD_A (local) DESCA[ LLD_ ] The leading dimension of the local
96 * array storing the local blocks of
97 * the distributed matrix A,
98 * IF( Lc( 1, N_A ) > 0 )
99 * LLD_A >= MAX( 1, Lr( 1, M_A ) )
100 * ELSE
101 * LLD_A >= 1.
102 *
103 * Let K be the number of rows of a matrix A starting at the global in-
104 * dex IA,i.e, A( IA:IA+K-1, : ). Lr( IA, K ) denotes the number of rows
105 * that the process of row coordinate MYROW ( 0 <= MYROW < NPROW ) would
106 * receive if these K rows were distributed over NPROW processes. If K
107 * is the number of columns of a matrix A starting at the global index
108 * JA, i.e, A( :, JA:JA+K-1, : ), Lc( JA, K ) denotes the number of co-
109 * lumns that the process MYCOL ( 0 <= MYCOL < NPCOL ) would receive if
110 * these K columns were distributed over NPCOL processes.
111 *
112 * The values of Lr() and Lc() may be determined via a call to the func-
113 * tion PB_Cnumroc:
114 * Lr( IA, K ) = PB_Cnumroc( K, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW )
115 * Lc( JA, K ) = PB_Cnumroc( K, JA, INB_A, NB_A, MYCOL, CSRC_A, NPCOL )
116 *
117 * Arguments
118 * =========
119 *
120 * TYPE (local input) pointer to a PBTYP_T structure
121 * On entry, TYPE is a pointer to a structure of type PBTYP_T,
122 * that contains type information (See pblas.h).
123 *
124 * CONJUG (global input) pointer to CHAR
125 * On entry, CONJUG specifies if this routine should return
126 * the conjugate subvector as follows:
127 * = 'N' or 'n': The initial subvector is returned,
128 * = 'Z' or 'z': The conjugate subvector is returned.
129 *
130 * ROWCOL (global input) pointer to CHAR
131 * On entry, ROWCOL specifies if this routine should return a
132 * row or column subvector replicated over the underlying subma-
133 * trix as follows:
134 * = 'R' or 'r': A row subvector is returned,
135 * = 'C' or 'c': A column subvector is returned.
136 *
137 * M (global input) INTEGER
138 * On entry, M specifies the number of rows of the underlying
139 * submatrix described by DESCA. M must be at least zero.
140 *
141 * N (global input) INTEGER
142 * On entry, N specifies the number of columns of the underlying
143 * submatrix described by DESCA. N must be at least zero.
144 *
145 * DESCA (global and local input) INTEGER array
146 * On entry, DESCA is an integer array of dimension DLEN_. This
147 * is the array descriptor for the matrix A.
148 *
149 * K (global input) INTEGER
150 * On entry, K specifies the length of the non-distributed di-
151 * mension of the subvector sub( X ). K must be at least zero.
152 *
153 * X (local input) pointer to CHAR
154 * On entry, X is an array of dimension (LLD_X, Kx), where LLD_X
155 * is at least MAX( 1, Lr( K, IX ) ) when XROC is 'R' or 'r'
156 * and MAX( 1, Lr( 1, IX+Lx-1 ) ) otherwise, and, Kx is at least
157 * Lc( 1, JX+Lx-1 ) when INCX = M_X and Lc( K, JX ) otherwise.
158 * Lx is N when ROWCOL = 'R' or 'r' and M otherwise. Before en-
159 * try, this array contains the local entries of the matrix X.
160 *
161 * IX (global input) INTEGER
162 * On entry, IX specifies X's global row index, which points to
163 * the beginning of the submatrix sub( X ).
164 *
165 * JX (global input) INTEGER
166 * On entry, JX specifies X's global column index, which points
167 * to the beginning of the submatrix sub( X ).
168 *
169 * DESCX (global and local input) INTEGER array
170 * On entry, DESCX is an integer array of dimension DLEN_. This
171 * is the array descriptor for the matrix X.
172 *
173 * XROC (global input) pointer to CHAR
174 * On entry, XROC specifies the orientation of the subvector
175 * sub( X ). When XROC is 'R' or 'r', sub( X ) is a row vector,
176 * and a column vector otherwise.
177 *
178 * XAPTR (local output) pointer to pointer to CHAR
179 * On exit, * XAPTR is an array containing the same data as the
180 * subvector sub( X ) which is replicated over the rows or co-
181 * lumns of the underlying matrix as specified by ROWCOL and
182 * DESCA.
183 *
184 * DXA (global and local output) INTEGER array
185 * On exit, DXA is a descriptor array of dimension DLEN_ descri-
186 * bing the data layout of the data pointed to by * XAPTR.
187 *
188 * XAFREE (local output) INTEGER
189 * On exit, XAFREE specifies if it has been possible to reuse
190 * the subvector sub( X ), i.e., if some dynamic memory was al-
191 * located for the data pointed to by *XAPTR or not. When XAFREE
192 * is zero, no dynamic memory was allocated. Otherwise, some dy-
193 * namic memory was allocated by this function that one MUST re-
194 * lease as soon as possible.
195 *
196 * -- Written on April 1, 1998 by
197 * Antoine Petitet, University of Tennessee, Knoxville 37996, USA.
198 *
199 * ---------------------------------------------------------------------
200 */
201 /*
202 * .. Local Scalars ..
203 */
204 char * top;
205 int AColSpan, ARowSpan, Acol, Aimb, Ainb, AisD, Amb, Amp, Anb,
206 Anq, Arow, Xcol, Xii, Ximb, Ximb1, Xinb, Xinb1, XisD, XisR,
207 XisRow, Xjj, Xld=1, Xmb, Xmp, Xnb, Xnq, Xrow, ctxt, mycol,
208 myrow, npcol, nprow;
209 /* ..
210 * .. Executable Statements ..
211 *
212 */
213 /*
214 * Initialize the output parameters to a default value
215 */
216 *XAFREE = 0;
217 *XAPTR = NULL;
218 /*
219 * Quick return if possible
220 */
221 if( ( M <= 0 ) || ( N <= 0 ) || ( K <= 0 ) )
222 {
223 if( Mupcase( ROWCOL[0] ) == CROW )
224 {
225 PB_Cdescset( DXA, K, N, 1, DESCA[INB_], 1, DESCA[NB_], DESCA[RSRC_],
226 DESCA[CSRC_], DESCA[CTXT_], 1 );
227 }
228 else
229 {
230 PB_Cdescset( DXA, M, K, DESCA[IMB_], 1, DESCA[MB_], 1, DESCA[RSRC_],
231 DESCA[CSRC_], DESCA[CTXT_], DESCA[LLD_] );
232 }
233 return;
234 }
235 /*
236 * Retrieve process grid information
237 */
238 Cblacs_gridinfo( ( ctxt = DESCX[CTXT_] ), &nprow, &npcol, &myrow, &mycol );
239 /*
240 * Retrieve sub( X )'s local information: Xii, Xjj, Xrow, Ycol
241 */
242 Minfog2l( IX, JX, DESCX, nprow, npcol, myrow, mycol, Xii, Xjj, Xrow, Xcol );
243 /*
244 * Is sub( X ) distributed or not, replicated or not ?
245 */
246 if( ( XisRow = ( Mupcase( XROC[0] ) == CROW ) ) != 0 )
247 {
248 XisD = ( ( Xcol >= 0 ) && ( npcol > 1 ) );
249 XisR = ( ( Xrow == -1 ) || ( nprow == 1 ) );
250 }
251 else
252 {
253 XisD = ( ( Xrow >= 0 ) && ( nprow > 1 ) );
254 XisR = ( ( Xcol == -1 ) || ( npcol == 1 ) );
255 }
256
257 Arow = DESCA[ RSRC_ ]; Acol = DESCA[ CSRC_ ];
258
259 if( Mupcase( ROWCOL[0] ) == CROW )
260 {
261 /*
262 * Want a row vector
263 */
264 Ainb = DESCA[ INB_ ]; Anb = DESCA[ NB_ ];
265 Mnumroc( Anq, N, 0, Ainb, Anb, mycol, Acol, npcol );
266 /*
267 * Does A spans multiples process rows ? It does if Arow < 0.
268 */
269 ARowSpan = ( Arow < 0 ) ||
270 Mspan( M, 0, DESCA[IMB_], DESCA[MB_], Arow, nprow );
271
272 if( XisRow && ( Mupcase( CONJUG[0] ) == CNOCONJG ) )
273 {
274 /*
275 * It is possible to reuse sub( X ) iff sub( X ) is already a row vector and
276 * the data does not need to be conjugated.
277 */
278 AisD = ( ( Acol >= 0 ) && ( npcol > 1 ) );
279
280 Xinb = DESCX[ INB_ ]; Xnb = DESCX[ NB_ ];
281 Mfirstnb( Xinb1, N, JX, Xinb, Xnb );
282 /*
283 * sub( X ) is aligned with A (reuse condition) iff both operands are not
284 * distributed, or both of them are distributed and start in the same process
285 * column and either N is smaller than the first blocksize of sub( X ) and A,
286 * or their column blocking factors match.
287 */
288 if( ( !AisD && !XisD ) ||
289 ( ( AisD && XisD ) &&
290 ( ( Acol == Xcol ) &&
291 ( ( ( Ainb >= N ) && ( Xinb1 >= N ) ) ||
292 ( ( Ainb == Xinb1 ) && ( Anb == Xnb ) ) ) ) ) )
293 {
294 /*
295 * sub( X ) is aligned with A
296 */
297 Ximb = DESCX[ IMB_ ]; Xmb = DESCX[ MB_ ];
298 Mfirstnb( Ximb1, K, IX, Ximb, Xmb );
299
300 if( XisR || ( !ARowSpan && ( Arow == Xrow ) ) )
301 {
302 /*
303 * If sub( X ) is replicated, or, A spans only one process row and either
304 * sub( X ) is replicated or resides in the same process row than A, then
305 * sub( X ) is already at the correct place.
306 */
307 if( Anq > 0 )
308 {
309 Xld = DESCX[ LLD_ ];
310 if( ARowSpan || ( myrow == Arow ) )
311 *XAPTR = Mptr( X, Xii, Xjj, Xld, TYPE->size );
312 }
313 else
314 {
315 Xld = 1;
316 }
317 MDescSet( DXA, K, N, K, Xinb1, 1, Xnb, ( ARowSpan ? -1 : Arow ),
318 Xcol, ctxt, Xld );
319 }
320 else if( ARowSpan )
321 {
322 /*
323 * Otherwise, we know that sub( X ) cannot be replicated, let suppose in
324 * addition that A spans all process rows. sub( X ) need simply to be broadcast
325 * over A.
326 */
327 if( myrow == Xrow )
328 {
329 Xld = DESCX[ LLD_ ];
330 if( Anq > 0 )
331 {
332 *XAPTR = Mptr( X, Xii, Xjj, Xld, TYPE->size );
333 top = PB_Ctop( &ctxt, BCAST, COLUMN, TOP_GET );
334 TYPE->Cgebs2d( ctxt, COLUMN, top, K, Anq, *XAPTR, Xld );
335 }
336 }
337 else
338 {
339 Xld = MAX( 1, K );
340 if( Anq > 0 )
341 {
342 *XAPTR = PB_Cmalloc( K * Anq * TYPE->size );
343 *XAFREE = 1;
344 top = PB_Ctop( &ctxt, BCAST, COLUMN, TOP_GET );
345 TYPE->Cgebr2d( ctxt, COLUMN, top, K, Anq, *XAPTR, Xld,
346 Xrow, mycol );
347 }
348 }
349 PB_Cdescset( DXA, K, N, K, Xinb1, 1, Xnb, -1, Xcol, ctxt, Xld );
350 }
351 else
352 {
353 /*
354 * Finally, sub( X ) is not replicated and A spans only one process row. There
355 * is no need to broadcast, a send/recv is sufficient.
356 */
357 if( myrow == Xrow )
358 {
359 Xld = DESCX[ LLD_ ];
360 if( Anq > 0 )
361 {
362 *XAPTR = Mptr( X, Xii, Xjj, Xld, TYPE->size );
363 TYPE->Cgesd2d( ctxt, K, Anq, *XAPTR, Xld, Arow, mycol );
364 }
365 }
366 else if( myrow == Arow )
367 {
368 Xld = MAX( 1, K );
369 if( Anq > 0 )
370 {
371 *XAPTR = PB_Cmalloc( K * Anq * TYPE->size );
372 *XAFREE = 1;
373 TYPE->Cgerv2d( ctxt, K, Anq, *XAPTR, Xld, Xrow, mycol );
374 }
375 }
376 PB_Cdescset( DXA, K, N, K, Xinb1, 1, Xnb, Arow, Xcol, ctxt,
377 Xld );
378 }
379 return;
380 }
381 }
382 /*
383 * sub( X ) cannot be reused, too bad ... redistribute
384 */
385 PB_Cdescset( DXA, K, N, K, Ainb, 1, Anb, ( ARowSpan ? -1 : Arow ), Acol,
386 ctxt, K );
387 Xmp = ( ARowSpan ? K : ( ( myrow == Arow ) ? K : 0 ) );
388 if( Xmp > 0 && Anq > 0 )
389 {
390 *XAPTR = PB_Cmalloc( Anq * Xmp * TYPE->size );
391 *XAFREE = 1;
392 }
393 if( XisRow )
394 {
395 PB_Cpaxpby( TYPE, CONJUG, K, N, TYPE->one, X, IX, JX, DESCX, XROC,
396 TYPE->zero, *XAPTR, 0, 0, DXA, ROW );
397 }
398 else
399 {
400 PB_Cpaxpby( TYPE, CONJUG, N, K, TYPE->one, X, IX, JX, DESCX, XROC,
401 TYPE->zero, *XAPTR, 0, 0, DXA, ROW );
402 }
403 }
404 else
405 {
406 /*
407 * Want a column vector
408 */
409 Aimb = DESCA[IMB_]; Amb = DESCA[MB_];
410 Mnumroc( Amp, M, 0, Aimb, Amb, myrow, Arow, nprow );
411 /*
412 * Does A spans multiples process columns ? It does if Acol < 0.
413 */
414 AColSpan = ( Acol < 0 ) ||
415 Mspan( N, 0, DESCA[INB_], DESCA[NB_], Acol, npcol );
416
417 if( !( XisRow ) && ( Mupcase( CONJUG[0] ) == CNOCONJG ) )
418 {
419 /*
420 * It is possible to reuse sub( X ) iff sub( X ) is already a column vector and
421 * the data does not need to be conjugated
422 */
423 AisD = ( ( Arow >= 0 ) && ( nprow > 1 ) );
424
425 Ximb = DESCX[ IMB_ ]; Xmb = DESCX[ MB_ ];
426 Mfirstnb( Ximb1, M, IX, Ximb, Xmb );
427 /*
428 * sub( X ) is aligned with A (reuse condition) iff both operands are not
429 * distributed, or both of them are distributed and start in the same process
430 * row and either M is smaller than the first blocksize of sub( X ) and A, or
431 * their row blocking factors match.
432 */
433 if( ( !AisD && !XisD ) ||
434 ( ( AisD && XisD ) &&
435 ( ( Arow == Xrow ) &&
436 ( ( ( Aimb >= M ) && ( Ximb1 >= M ) ) ||
437 ( ( Aimb == Ximb1 ) && ( Amb == Xmb ) ) ) ) ) )
438 {
439 /*
440 * sub( X ) is aligned with A
441 */
442 Xinb = DESCX[ INB_ ]; Xnb = DESCX[ NB_ ];
443 Mfirstnb( Xinb1, K, JX, Xinb, Xnb );
444
445 if( XisR || ( !AColSpan && ( Acol == Xcol ) ) )
446 {
447 /*
448 * If sub( X ) is replicated, or, A spans only one process column and either
449 * sub( X ) is replicated or resides in the same process columns than A, then
450 * sub( X ) is already at the correct place.
451 */
452 if( Amp > 0 )
453 {
454 Xld = DESCX[ LLD_ ];
455 if( AColSpan || ( mycol == Acol ) )
456 *XAPTR = Mptr( X, Xii, Xjj, Xld, TYPE->size );
457 }
458 else
459 {
460 Xld = 1;
461 }
462 MDescSet( DXA, M, K, Ximb1, K, Xmb, 1, Xrow,
463 ( AColSpan ? -1 : Acol ), ctxt, Xld );
464 }
465 else if( AColSpan )
466 {
467 /*
468 * Otherwise, we know that sub( X ) is not be replicated, let suppose in
469 * addition that A spans all process columns. sub( X ) need simply to be
470 * broadcast over A.
471 */
472 if( mycol == Xcol )
473 {
474 Xld = DESCX[ LLD_ ];
475 if( Amp > 0 )
476 {
477 *XAPTR = Mptr( X, Xii, Xjj, Xld, TYPE->size );
478 top = PB_Ctop( &ctxt, BCAST, ROW, TOP_GET );
479 TYPE->Cgebs2d( ctxt, ROW, top, Amp, K, *XAPTR, Xld );
480 }
481 }
482 else
483 {
484 Xld = MAX( 1, Amp );
485 if( Amp > 0 )
486 {
487 *XAPTR = PB_Cmalloc( Amp * K * TYPE->size );
488 *XAFREE = 1;
489 top = PB_Ctop( &ctxt, BCAST, ROW, TOP_GET );
490 TYPE->Cgebr2d( ctxt, ROW, top, Amp, K, *XAPTR, Xld, myrow,
491 Xcol );
492 }
493 }
494 PB_Cdescset( DXA, M, K, Ximb1, K, Xmb, 1, Xrow, -1, ctxt, Xld );
495 }
496 else
497 {
498 /*
499 * Finally, sub( X ) is not replicated and A spans only one process column.
500 * There is no need to broadcast, a send/recv is sufficient.
501 */
502 if( mycol == Xcol )
503 {
504 Xld = DESCX[ LLD_ ];
505 if( Amp > 0 )
506 {
507 *XAPTR = Mptr( X, Xii, Xjj, Xld, TYPE->size );
508 TYPE->Cgesd2d( ctxt, Amp, K, *XAPTR, Xld, myrow, Acol );
509 }
510 }
511 else if( mycol == Acol )
512 {
513 Xld = MAX( 1, Amp );
514 if( Amp > 0 )
515 {
516 *XAPTR = PB_Cmalloc( Amp * K * TYPE->size );
517 *XAFREE = 1;
518 TYPE->Cgerv2d( ctxt, Amp, K, *XAPTR, Xld, myrow, Xcol );
519 }
520 }
521 PB_Cdescset( DXA, M, K, Ximb1, K, Xmb, 1, Xrow, Acol, ctxt,
522 Xld );
523 }
524 return;
525 }
526 }
527 /*
528 * sub( X ) cannot be reused, too bad ... redistribute
529 */
530 PB_Cdescset( DXA, M, K, Aimb, K, Amb, 1, Arow, ( AColSpan ? -1 : Acol ),
531 ctxt, MAX( 1, Amp ) );
532 Xnq = ( AColSpan ? K : ( ( mycol == Acol ) ? K : 0 ) );
533 if( Xnq > 0 && Amp > 0 )
534 {
535 *XAPTR = PB_Cmalloc( Amp * Xnq * TYPE->size );
536 *XAFREE = 1;
537 }
538 if( XisRow )
539 {
540 PB_Cpaxpby( TYPE, CONJUG, K, M, TYPE->one, X, IX, JX, DESCX, XROC,
541 TYPE->zero, *XAPTR, 0, 0, DXA, COLUMN );
542 }
543 else
544 {
545 PB_Cpaxpby( TYPE, CONJUG, M, K, TYPE->one, X, IX, JX, DESCX, XROC,
546 TYPE->zero, *XAPTR, 0, 0, DXA, COLUMN );
547 }
548 }
549 /*
550 * End of PB_CInV
551 */
552 }
553