SpM Handbook 1.2.4
Loading...
Searching...
No Matches
spm_gather.c
Go to the documentation of this file.
1/**
2 *
3 * @file spm_gather.c
4 *
5 * SParse Matrix gather routine.
6 *
7 * @copyright 2020-2024 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
8 * Univ. Bordeaux. All rights reserved.
9 *
10 * @version 1.2.4
11 * @author Tony Delarue
12 * @author Mathieu Faverge
13 * @author Alycia Lisito
14 * @date 2024-06-27
15 *
16 **/
17#include "common.h"
18
19#if !defined(SPM_WITH_MPI)
20#error "This file should not be compiled if MPI support is not enabled (SPM_WITH_MPI)"
21#endif
22
23/**
24 *******************************************************************************
25 *
26 * @brief Initialize the gather structures
27 *
28 *******************************************************************************
29 *
30 * @param[in] spm
31 * The spm to gather.
32 *
33 *******************************************************************************
34 *
35 * @return The array of triplets { n, nnz, nnzexp } for all nodes.
36 *
37 *******************************************************************************/
38static inline int *
40{
41 int *allcounts = NULL;
42 int values[3] = { spm->n, spm->nnz, spm->nnzexp };
43
44 allcounts = malloc( 3 * spm->clustnbr * sizeof(int) );
45
46 /* Collect the number of elements per node */
47 MPI_Allgather( values, 3, MPI_INT,
48 allcounts, 3, MPI_INT, spm->comm );
49
50 return allcounts;
51}
52
53/**
54 *******************************************************************************
55 *
56 * @brief Check if a matrix has been scattered continuously or not.
57 *
58 *******************************************************************************
59 *
60 * @param[in] spm
61 * The spm to gather.
62 *
63 * @param[in] allcounts
64 * The array of triplets { n, nnz, nnzexp } for all nodes.
65 *
66 *******************************************************************************
67 *
68 * @retval 0 if the matrix is not scattered continuously
69 * @retval !0 if the matrix is scattered continuously
70 *
71 *******************************************************************************/
72static inline int
74 const int *allcounts )
75{
76 spm_int_t baseval = spm->baseval;
77 spm_int_t shift = baseval;
78 spm_int_t i;
79 int rc;
80 int error = 0;
81
82 assert( spm->fmttype != SpmIJV );
83
84 for ( i = 0; i < spm->clustnum; i++ ) {
85 shift += allcounts[i * 3];
86 }
87
88 /* Check its local loc2glob */
89 for ( i = 0; i < spm->n; i++ ) {
90 if ( spm->loc2glob[i] != (shift + i) ) {
91 error = 1;
92 break;
93 }
94 }
95
96 MPI_Allreduce( &error, &rc, 1, MPI_INT,
97 MPI_SUM, spm->comm );
98
99 return !rc;
100}
101
102/**
103 *******************************************************************************
104 *
105 * @brief Update a gathered compressed array to shift the indices accordingly to
106 * the distribution.
107 *
108 *******************************************************************************
109 *
110 * @param[in] spm
111 * The original scattered spm
112 *
113 * @param[in] colptr
114 * The pointer to the gathered compressed array (spm->colptr, or
115 * spm->rowptr) of the new gathered spm.
116 *
117 * @param[in] recvdispls
118 * The array of reception displacements for the n values.
119 *
120 * @param[in] recvcounts
121 * The array of reception count in terms of nnz.
122 *
123 *******************************************************************************/
124static inline void
126 spm_int_t *colptr,
127 int *recvdispls,
128 int *recvcounts )
129{
130 int c;
131 spm_int_t to_add = 0;
132 /*
133 * We need to update the compressed array to match the gathered array
134 *
135 * rowptr[ 0:N_0-1] += 0
136 * rowptr[N_0:N_1-1] += NNZ_0
137 * rowptr[N_1:N_2-1] += NNZ_0 + NNZ_1
138 * ...
139 * rowptr[N] += NNZ
140 */
141 for ( c=1; c<spm->clustnbr; c++ ) {
142 /* Let's start shifting the value after the first array */
143 spm_int_t shift = recvdispls[c];
144 spm_int_t end = ( c < (spm->clustnbr - 1) ) ? recvdispls[c+1] : spm->gN;
145 spm_int_t i;
146
147 to_add += recvcounts[c-1];
148 for ( i=shift; i<end; i++ ) {
149 colptr[i] += to_add;
150 }
151 }
152 assert( to_add + recvcounts[spm->clustnbr-1] == spm->gnnz );
153 /* Set the last value */
154 colptr[spm->gN] = colptr[0] + spm->gnnz;
155}
156
157/**
158 *******************************************************************************
159 *
160 * @brief Gather a distributed Sparse Matrix on the root node(s) in CSC/CSR format.
161 *
162 *******************************************************************************
163 *
164 * @param[in] oldspm
165 * The distributed sparse matrix to gather.
166 *
167 * @param[in] newspm
168 * The new gathered spm. NULL if not root.
169 *
170 * @param[in] root
171 * The root node that gather the final sparse matrix.
172 * If -1, all nodes gather a copy of the matrix.
173 *
174 * @param[in] allcounts
175 * The array of the triplets {n, nnz, nnzexp}.
176 *
177 *******************************************************************************/
178static inline void
180 spmatrix_t *newspm,
181 int root,
182 int *allcounts )
183{
184 const spm_int_t *oldcol = ( oldspm->fmttype == SpmCSC ) ? oldspm->colptr : oldspm->rowptr;
185 const spm_int_t *oldrow = ( oldspm->fmttype == SpmCSC ) ? oldspm->rowptr : oldspm->colptr;
186 const char *oldval = oldspm->values;
187 spm_int_t *newcol = NULL;
188 spm_int_t *newrow = NULL;
189 char *newval = NULL;
190
191 int *recvcounts = NULL;
192 int *recvdispls = NULL;
193 int recv = ( root == -1 ) || ( root == oldspm->clustnum );
194 int c;
195
196 assert( ((newspm != NULL) && recv) ||
197 ((newspm == NULL) && !recv) );
198
199 /*
200 * First step, let's gather the compressed array
201 */
202 {
203 int n = oldspm->n;
204 if ( recv ) {
205 recvcounts = malloc( oldspm->clustnbr * sizeof(int) );
206 recvdispls = malloc( oldspm->clustnbr * sizeof(int) );
207
208 recvdispls[0] = 0;
209 recvcounts[0] = allcounts[0]; /* n */
210 for( c=1; c<oldspm->clustnbr; c++ ) {
211 recvdispls[c] = recvdispls[c-1] + recvcounts[c-1];
212 recvcounts[c] = allcounts[3 * c];
213 }
214
215 /* Initialize the new pointers */
216 assert( newspm );
217 newcol = ( newspm->fmttype == SpmCSC ) ? newspm->colptr : newspm->rowptr;
218 newrow = ( newspm->fmttype == SpmCSC ) ? newspm->rowptr : newspm->colptr;
219 newval = newspm->values;
220 }
221
222 /* Gather the colptr */
223 if ( root == -1 ) {
224 MPI_Allgatherv( oldcol, n, SPM_MPI_INT,
225 newcol, recvcounts, recvdispls, SPM_MPI_INT, oldspm->comm );
226 }
227 else {
228 MPI_Gatherv( oldcol, n, SPM_MPI_INT,
229 newcol, recvcounts, recvdispls, SPM_MPI_INT, root, oldspm->comm );
230 }
231
232 if ( recv ) {
233 /* recvdispls : n, recvcnt : nnz */
234 recvcounts[0] = allcounts[1]; /* nnz */
235 for( c=1; c<oldspm->clustnbr; c++ ) {
236 recvcounts[c] = allcounts[3 * c + 1];
237 }
238 spm_gather_csx_update( oldspm, newcol, recvdispls, recvcounts );
239 }
240 }
241
242 /*
243 * Second step, let's gather the non compressed array
244 */
245 {
246 int nnz = oldspm->nnz;
247
248 if ( recv ) {
249 recvdispls[0] = 0;
250 recvcounts[0] = allcounts[1]; /* nnz */
251 for( c=1; c<oldspm->clustnbr; c++ ) {
252 recvdispls[c] = recvdispls[c-1] + recvcounts[c-1];
253 recvcounts[c] = allcounts[3 * c + 1];
254 }
255 }
256
257 /* Gather the arrays */
258 if ( root == -1 ) {
259 MPI_Allgatherv( oldrow, nnz, SPM_MPI_INT,
260 newrow, recvcounts, recvdispls, SPM_MPI_INT, oldspm->comm );
261 }
262 else {
263 MPI_Gatherv( oldrow, nnz, SPM_MPI_INT,
264 newrow, recvcounts, recvdispls, SPM_MPI_INT,
265 root, oldspm->comm );
266 }
267 }
268
269 /*
270 * Third step, let's gather the values array
271 */
272 if ( oldspm->flttype != SpmPattern ) {
273 MPI_Datatype valtype = spm_get_datatype( oldspm );
274 int nnzexp = oldspm->nnzexp;
275
276 if ( recv ) {
277 recvdispls[0] = 0;
278 recvcounts[0] = allcounts[2]; /* nnzexp */
279 for( c=1; c<oldspm->clustnbr; c++ ) {
280 recvdispls[c] = recvdispls[c-1] + recvcounts[c-1];
281 recvcounts[c] = allcounts[3 * c + 2];
282 }
283 }
284
285 if ( root == -1 ) {
286 MPI_Allgatherv( oldval, nnzexp, valtype,
287 newval, recvcounts, recvdispls, valtype, oldspm->comm );
288 }
289 else {
290 MPI_Gatherv( oldval, nnzexp, valtype,
291 newval, recvcounts, recvdispls, valtype,
292 root, oldspm->comm );
293 }
294 }
295
296 /* Cleanup */
297 if ( recv ) {
298 free( recvcounts );
299 free( recvdispls );
300 }
301}
302
303/**
304 *******************************************************************************
305 *
306 * @brief Gather a distributed Sparse Matrix on the root node(s) in IJV format.
307 *
308 *******************************************************************************
309 *
310 * @param[in] oldspm
311 * The distributed sparse matrix to gather.
312 *
313 * @param[in] newspm
314 * The new gathered spm. NULL if not root.
315 *
316 * @param[in] root
317 * The root node that gather the final sparse matrix.
318 * If -1, all nodes gather a copy of the matrix.
319 *
320 * @param[in] allcounts
321 * The array of the triplets {n, nnz, nnzexp}.
322 *
323 *******************************************************************************/
324static inline void
326 spmatrix_t *newspm,
327 int root,
328 const int *allcounts )
329{
330 MPI_Datatype valtype = spm_get_datatype( oldspm );
331 int *recvcounts = NULL;
332 int *recvdispls = NULL;
333 int nnz = oldspm->nnz;
334 int nnzexp = oldspm->nnzexp;
335 int recv = ( root == -1 ) || ( root == oldspm->clustnum );
336
337 assert( ((newspm != NULL) && recv) ||
338 ((newspm == NULL) && !recv) );
339
340 if ( recv ) {
341 int c;
342 recvcounts = malloc( oldspm->clustnbr * sizeof(int) );
343 recvdispls = malloc( oldspm->clustnbr * sizeof(int) );
344
345 recvdispls[0] = 0;
346 recvcounts[0] = allcounts[1];
347 for( c=1; c<oldspm->clustnbr; c++ ) {
348 recvdispls[c] = recvdispls[c-1] + recvcounts[c-1];
349 recvcounts[c] = allcounts[3 * c + 1];
350 }
351 }
352
353 /* Gather the indices arrays */
354 if ( root == -1 ) {
355 MPI_Allgatherv( oldspm->colptr, nnz, SPM_MPI_INT,
356 newspm->colptr, recvcounts, recvdispls, SPM_MPI_INT, oldspm->comm );
357 MPI_Allgatherv( oldspm->rowptr, nnz, SPM_MPI_INT,
358 newspm->rowptr, recvcounts, recvdispls, SPM_MPI_INT, oldspm->comm );
359 }
360 else {
361 spm_int_t *newcol = recv ? newspm->colptr : NULL;
362 spm_int_t *newrow = recv ? newspm->rowptr : NULL;
363
364 MPI_Gatherv( oldspm->colptr, nnz, SPM_MPI_INT,
365 newcol, recvcounts, recvdispls, SPM_MPI_INT,
366 root, oldspm->comm );
367
368 MPI_Gatherv( oldspm->rowptr, nnz, SPM_MPI_INT,
369 newrow, recvcounts, recvdispls, SPM_MPI_INT,
370 root, oldspm->comm );
371 }
372
373 /* Gather the values */
374 if ( oldspm->flttype != SpmPattern )
375 {
376 /* Update recvcounts and recvdispls arrays if needed to use nnzexp */
377 if ( recv && ( oldspm->dof != 1 ) ) {
378 int c;
379
380 recvdispls[0] = 0;
381 recvcounts[0] = allcounts[2];
382 for( c=1; c<oldspm->clustnbr; c++ ) {
383 recvdispls[c] = recvdispls[c-1] + recvcounts[c-1];
384 recvcounts[c] = allcounts[3 * c + 2];
385 }
386 }
387
388 if ( root == -1 ) {
389 MPI_Allgatherv( oldspm->values, nnzexp, valtype,
390 newspm->values, recvcounts, recvdispls, valtype, oldspm->comm );
391 }
392 else {
393 char *newval = recv ? newspm->values : NULL;
394 MPI_Gatherv( oldspm->values, nnzexp, valtype,
395 newval, recvcounts, recvdispls, valtype,
396 root, oldspm->comm );
397 }
398 }
399
400 if ( recv ) {
401 free( recvcounts );
402 free( recvdispls );
403 }
404}
405
406/**
407 *******************************************************************************
408 *
409 * @brief Gather a distributed Sparse Matrix on a single node.
410 *
411 * The matrix must be distributed continously on each node to be gathered. If
412 * it's not the case, please apply first a permutation where the loc2glob array
413 * corresponds to the inverse permutation.
414 *
415 *******************************************************************************
416 *
417 * @param[in] oldspm
418 * The distributed sparse matrix to gather.
419 *
420 * @param[in] root
421 * The root node that gather the final sparse matrix.
422 * If -1, all nodes gather a copy of the matrix.
423 *
424 * @param[inout] newspm
425 * On entry, the allocated spm structure if we are on one root of the
426 * gather. Not referenced otherwise.
427 * On exit, contains the gathered spm on root node(s).
428 *
429 *******************************************************************************
430 *
431 * @retval SPM_SUCCESS on success.
432 * @retval SPM_ERR_BADPARAMETER on incorrect parameter.
433 *
434 *******************************************************************************/
435int
436spmGather( const spmatrix_t *oldspm,
437 int root,
438 spmatrix_t *newspm )
439{
440 spmatrix_t *spmd = NULL;
441 int *allcounts;
442
443 assert( oldspm->replicated != -1 );
444 if ( ( oldspm->clustnbr == 1 ) || ( oldspm->replicated ) ) {
445 if ( newspm == NULL ) {
446 spm_print_error("spmGather: Incorrect call with newspm == NULL\n");
448 }
449 spmCopy( oldspm, newspm );
450 newspm->replicated = 1;
451 return SPM_SUCCESS;
452 }
453
454 assert( !oldspm->replicated );
455 assert( oldspm->loc2glob != NULL );
456 assert( ( root >= -1 ) && ( root < oldspm->clustnbr ) );
457
458 allcounts = spm_gather_init( oldspm );
459
460 if ( oldspm->fmttype != SpmIJV ) {
461 /* Check if the matrix is distributed continuously */
462 int continuous = spm_gather_check( oldspm, allcounts );
463
464 if ( continuous ) {
465 spmd = (spmatrix_t *)oldspm;
466 }
467 else {
468 spm_int_t *newl2g, new_n;
469 spmd = malloc( sizeof(spmatrix_t) );
470
471 /* Redistribute the matrix continuously */
472 new_n = spm_create_loc2glob_continuous( oldspm, &newl2g );
473 spmRedistribute( oldspm, new_n, newl2g, spmd );
474 free( newl2g );
475 free( allcounts );
476
477 /* Reinitialize allcounts */
478 allcounts = spm_gather_init( spmd );
479 assert( spm_gather_check(spmd, allcounts) );
480 }
481 }
482 else {
483 spmd = (spmatrix_t *)oldspm;
484 }
485 assert( spmd != NULL );
486
487 if ( (root == -1) ||
488 (root == oldspm->clustnum) )
489 {
490 if ( newspm == NULL ) {
491 spm_print_error("spmGather: Incorrect call with newspm == NULL\n");
493 }
494 spmInit( newspm );
495 memcpy( newspm, spmd, sizeof(spmatrix_t) );
496
497 /* Compute specific informations */
498 newspm->baseval = spmd->baseval;
499 newspm->n = spmd->gN;
500 newspm->nexp = spmd->gNexp;
501 newspm->nnz = spmd->gnnz;
502 newspm->nnzexp = spmd->gnnzexp;
503
504 newspm->replicated = 1;
505
506 newspm->dofs = NULL;
507 newspm->colptr = NULL;
508 newspm->rowptr = NULL;
509 newspm->values = NULL;
510
511 newspm->loc2glob = NULL;
512 newspm->glob2loc = NULL;
513
514 if ( root != -1 ) {
515 newspm->clustnbr = 1;
516 newspm->comm = MPI_COMM_SELF;
517 }
518
519 spmAlloc( newspm );
520 if ( newspm->dof < 1 ) {
521 memcpy( newspm->dofs, spmd->dofs, ((size_t)(newspm->gN + 1)) * sizeof(spm_int_t) );
522 }
523 }
524
525 switch( spmd->fmttype ) {
526 case SpmCSC:
527 case SpmCSR:
528 spm_gather_csx( spmd, newspm, root, allcounts );
529 break;
530 case SpmIJV:
531 default:
532 spm_gather_ijv( spmd, newspm, root, allcounts );
533 }
534
535 free( allcounts );
536
537 /* Free the possibly redistributed spm */
538 if ( spmd != oldspm ) {
539 spmExit(spmd);
540 free(spmd);
541 }
542
543 return SPM_SUCCESS;
544}
545
546/**
547 *******************************************************************************
548 *
549 * @brief This routine performs a allgather of a distributed spm in place.
550 *
551 *******************************************************************************
552 *
553 * @param[in] spm
554 * On entry, the distributed spm.
555 * On exit, the gathered (replicated) spm.
556 *
557 *******************************************************************************
558 *
559 * @retval 1 if the spm has been gathered, 0 if untouched
560 * Note that untouched is not necessarily a failure if the spm was not
561 * distributed.
562 *
563 ********************************************************************************/
564int
566{
567 spmatrix_t newspm;
568 int rc;
569
570 assert( spm->replicated != -1 );
571 if ( spm->replicated || (spm->clustnbr == 1) ) {
572 return 0;
573 }
574
575 rc = spmGather( spm, -1, &newspm );
576
577 if ( rc == SPM_SUCCESS ) {
578 spmExit( spm );
579 memcpy( spm, &newspm, sizeof( spmatrix_t ) );
580 return 1;
581 }
582 else {
583 return 0;
584 }
585}
@ SPM_ERR_BADPARAMETER
Definition const.h:89
@ SPM_SUCCESS
Definition const.h:82
@ SpmPattern
Definition const.h:62
@ SpmCSC
Definition const.h:73
@ SpmCSR
Definition const.h:74
@ SpmIJV
Definition const.h:75
spm_coeftype_t flttype
Definition spm.h:67
spm_int_t * dofs
Definition spm.h:85
spm_int_t nexp
Definition spm.h:78
void * values
Definition spm.h:92
int clustnbr
Definition spm.h:96
spm_int_t nnz
Definition spm.h:75
spm_int_t nnzexp
Definition spm.h:80
spm_int_t * rowptr
Definition spm.h:90
spm_int_t gnnz
Definition spm.h:74
spm_fmttype_t fmttype
Definition spm.h:69
spm_int_t dof
Definition spm.h:82
spm_int_t n
Definition spm.h:73
spm_int_t * loc2glob
Definition spm.h:91
SPM_Comm comm
Definition spm.h:97
spm_int_t baseval
Definition spm.h:71
spm_int_t * glob2loc
Definition spm.h:94
spm_int_t gnnzexp
Definition spm.h:79
spm_int_t * colptr
Definition spm.h:89
spm_int_t gNexp
Definition spm.h:77
int clustnum
Definition spm.h:95
spm_int_t gN
Definition spm.h:72
int replicated
Definition spm.h:98
int spmGather(const spmatrix_t *oldspm, int root, spmatrix_t *newspm)
Gather a distributed Sparse Matrix on a single node.
Definition spm_gather.c:436
void spmAlloc(spmatrix_t *spm)
Allocate the arrays of an spm structure.
Definition spm.c:175
void spmExit(spmatrix_t *spm)
Cleanup the spm structure but do not free the spm pointer.
Definition spm.c:224
#define SPM_MPI_INT
The MPI type associated to spm_int_t.
Definition datatypes.h:72
int spm_int_t
The main integer datatype used in spm arrays.
Definition datatypes.h:70
int spmRedistribute(const spmatrix_t *spm, spm_int_t new_n, const spm_int_t *newl2g, spmatrix_t *newspm)
Create a copy of a given SPM with a new distribution corresponding to the given loc2glob array.
int spmGatherInPlace(spmatrix_t *spm)
This routine performs a allgather of a distributed spm in place.
Definition spm_gather.c:565
void spmCopy(const spmatrix_t *spm_in, spmatrix_t *spm_out)
Create a copy of the spm.
Definition spm.c:969
void spmInit(spmatrix_t *spm)
Init the spm structure.
Definition spm.c:152
The sparse matrix data structure.
Definition spm.h:64
spm_int_t spm_create_loc2glob_continuous(const spmatrix_t *spm, spm_int_t **l2g_ptr)
Generate a continuous loc2glob array on each node.
Definition spm.c:1841
static void spm_gather_ijv(const spmatrix_t *oldspm, spmatrix_t *newspm, int root, const int *allcounts)
Gather a distributed Sparse Matrix on the root node(s) in IJV format.
Definition spm_gather.c:325
static int * spm_gather_init(const spmatrix_t *spm)
Initialize the gather structures.
Definition spm_gather.c:39
static void spm_gather_csx_update(const spmatrix_t *spm, spm_int_t *colptr, int *recvdispls, int *recvcounts)
Update a gathered compressed array to shift the indices accordingly to the distribution.
Definition spm_gather.c:125
static void spm_gather_csx(const spmatrix_t *oldspm, spmatrix_t *newspm, int root, int *allcounts)
Gather a distributed Sparse Matrix on the root node(s) in CSC/CSR format.
Definition spm_gather.c:179
static int spm_gather_check(const spmatrix_t *spm, const int *allcounts)
Check if a matrix has been scattered continuously or not.
Definition spm_gather.c:73