1 /*============================================================================
2 * Manipulation of global indexed list
3 *===========================================================================*/
4
5 /*
6 This file is part of Code_Saturne, a general-purpose CFD tool.
7
8 Copyright (C) 1998-2021 EDF S.A.
9
10 This program is free software; you can redistribute it and/or modify it under
11 the terms of the GNU General Public License as published by the Free Software
12 Foundation; either version 2 of the License, or (at your option) any later
13 version.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
18 details.
19
20 You should have received a copy of the GNU General Public License along with
21 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
22 Street, Fifth Floor, Boston, MA 02110-1301, USA.
23 */
24
25 /*----------------------------------------------------------------------------*/
26
27 #include "cs_defs.h"
28
29 /*----------------------------------------------------------------------------
30 * Standard C library headers
31 *---------------------------------------------------------------------------*/
32
33 #include <assert.h>
34 #include <stdio.h>
35 #include <string.h>
36
37 /*----------------------------------------------------------------------------
38 * Local headers
39 *---------------------------------------------------------------------------*/
40
41 #include "bft_mem.h"
42
43 #include "fvm_io_num.h"
44
45 #include "cs_all_to_all.h"
46 #include "cs_block_dist.h"
47 #include "cs_join_util.h"
48 #include "cs_order.h"
49 #include "cs_search.h"
50 #include "cs_sort.h"
51
52 /*----------------------------------------------------------------------------
53 * Header for the current file
54 *---------------------------------------------------------------------------*/
55
56 #include "cs_join_set.h"
57
58 /*---------------------------------------------------------------------------*/
59
60 BEGIN_C_DECLS
61
62 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
63
64 /*============================================================================
65 * Macro and type definitions
66 *===========================================================================*/
67
68 /*============================================================================
69 * Private function definitions
70 *===========================================================================*/
71
72 /*----------------------------------------------------------------------------
73 * Sort an array "a" and apply the sort to its associated array "b"
74 * (global numbering).
75 *
76 * Sort is realized using a shell sort (Knuth algorithm).
77 *
78 * parameters:
79 * l <-- left bound
80 * r <-- right bound
81 * a <-> array to sort
82 * b <-> associated array
83 *---------------------------------------------------------------------------*/
84
85 static void
_coupled_adapted_gnum_shellsort(cs_lnum_t l,cs_lnum_t r,cs_gnum_t a[],cs_gnum_t b[])86 _coupled_adapted_gnum_shellsort(cs_lnum_t l,
87 cs_lnum_t r,
88 cs_gnum_t a[],
89 cs_gnum_t b[])
90 {
91 cs_lnum_t i, start;
92 cs_gnum_t ref;
93
94 cs_lnum_t size = r - l;
95
96 if (size == 0)
97 return;
98
99 cs_sort_coupled_gnum_shell(l, r, a, b);
100
101 /* Order b array for each sub-list where a[] is constant */
102
103 start = l;
104 while (start < r) {
105
106 ref = a[start];
107 for (i = start; i < r && ref == a[i]; i++);
108
109 cs_sort_gnum_shell(start, i, b);
110 start = i;
111
112 }
113
114 }
115
116 /*----------------------------------------------------------------------------
117 * Test if an array of local numbers with stride 2 is lexicographically
118 * ordered.
119 *
120 * parameters:
121 * number <-- array of all entity numbers (number of entity i
122 * given by number[i] or number[list[i] - 1])
123 * nb_ent <-- number of entities considered
124 *
125 * returns:
126 * 1 if ordered, 0 otherwise.
127 *----------------------------------------------------------------------------*/
128
129 static int
_order_local_test_s2(const cs_lnum_t number[],size_t n_elts)130 _order_local_test_s2(const cs_lnum_t number[],
131 size_t n_elts)
132 {
133 size_t i = 0;
134
135 assert (number != NULL || n_elts == 0);
136
137 for (i = 1 ; i < n_elts ; i++) {
138 size_t i_prev, k;
139 bool unordered = false;
140 i_prev = i-1;
141 for (k = 0; k < 2; k++) {
142 if (number[i_prev*2 + k] < number[i*2 + k])
143 break;
144 else if (number[i_prev*2 + k] > number[i*2 + k])
145 unordered = true;
146 }
147 if (unordered == true)
148 break;
149 }
150
151 if (i == n_elts || n_elts == 0)
152 return 1;
153 else
154 return 0;
155 }
156
157 /*----------------------------------------------------------------------------
158 * Descend binary tree for the lexicographical ordering of an local numbering
159 * array of stride 2.
160 *
161 * parameters:
162 * number <-- pointer to numbers of entities that should be ordered.
163 * level <-- level of the binary tree to descend
164 * n_elts <-- number of entities in the binary tree to descend
165 * order <-> ordering array
166 *----------------------------------------------------------------------------*/
167
168 inline static void
_order_descend_tree_s2(const cs_lnum_t number[],size_t level,const size_t n_elts,cs_lnum_t order[])169 _order_descend_tree_s2(const cs_lnum_t number[],
170 size_t level,
171 const size_t n_elts,
172 cs_lnum_t order[])
173 {
174 size_t i_save, i1, i2, j, lv_cur;
175
176 i_save = (size_t)(order[level]);
177
178 while (level <= (n_elts/2)) {
179
180 lv_cur = (2*level) + 1;
181
182 if (lv_cur < n_elts - 1) {
183
184 i1 = (size_t)(order[lv_cur+1]);
185 i2 = (size_t)(order[lv_cur]);
186
187 for (j = 0; j < 2; j++) {
188 if (number[i1*2 + j] != number[i2*2 + j])
189 break;
190 }
191
192 if (j < 2) {
193 if (number[i1*2 + j] > number[i2*2 + j])
194 lv_cur++;
195 }
196
197 }
198
199 if (lv_cur >= n_elts) break;
200
201 i1 = i_save;
202 i2 = (size_t)(order[lv_cur]);
203
204 for (j = 0; j < 2; j++) {
205 if (number[i1*2 + j] != number[i2*2 + j])
206 break;
207 }
208
209 if (j == 2) break;
210 if (number[i1*2 + j] >= number[i2*2 + j]) break;
211
212 order[level] = order[lv_cur];
213 level = lv_cur;
214
215 }
216
217 order[level] = i_save;
218 }
219
220 /*----------------------------------------------------------------------------
221 * Order array of local numbers with stride 2 lexicographically.
222 *
223 * parameters:
224 * number <-- array of entity numbers
225 * order --> pre-allocated ordering table
226 * n_elts <-- number of entities considered
227 *----------------------------------------------------------------------------*/
228
229 static void
_order_local_s2(const cs_lnum_t number[],cs_lnum_t order[],const size_t n_elts)230 _order_local_s2(const cs_lnum_t number[],
231 cs_lnum_t order[],
232 const size_t n_elts)
233 {
234 size_t i;
235 cs_lnum_t o_save;
236
237 assert (number != NULL || n_elts == 0);
238
239 /* Initialize ordering array */
240
241 for (i = 0 ; i < n_elts ; i++)
242 order[i] = i;
243
244 if (n_elts < 2)
245 return;
246
247 /* Create binary tree */
248
249 i = (n_elts / 2) ;
250 do {
251 i--;
252 _order_descend_tree_s2(number, i, n_elts, order);
253 } while (i > 0);
254
255 /* Sort binary tree */
256
257 for (i = n_elts - 1 ; i > 0 ; i--) {
258 o_save = order[0];
259 order[0] = order[i];
260 order[i] = o_save;
261 _order_descend_tree_s2(number, 0, i, order);
262 }
263 }
264
265 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
266
267 /*============================================================================
268 * Public function definitions
269 *===========================================================================*/
270
271 /*----------------------------------------------------------------------------
272 * Allocate a resizable array.
273 *
274 * parameters:
275 * max_size <-- initial number of elements to allocate
276 *
277 * returns:
278 * pointer to a new alloacted resizable array
279 *---------------------------------------------------------------------------*/
280
281 cs_join_rset_t *
cs_join_rset_create(cs_lnum_t max_size)282 cs_join_rset_create(cs_lnum_t max_size)
283 {
284 cs_join_rset_t *new_set = NULL;
285
286 if (max_size > 0) {
287
288 BFT_MALLOC(new_set, 1, cs_join_rset_t);
289
290 new_set->n_max_elts = max_size;
291 new_set->n_elts = 0;
292
293 BFT_MALLOC(new_set->array, max_size, cs_lnum_t);
294
295 }
296
297 return new_set;
298 }
299
300 /*----------------------------------------------------------------------------
301 * Destroy a cs_join_rset_t structure.
302 *
303 * parameters:
304 * set <-- pointer to pointer to the cs_join_rset_t structure to destroy
305 *---------------------------------------------------------------------------*/
306
307 void
cs_join_rset_destroy(cs_join_rset_t ** set)308 cs_join_rset_destroy(cs_join_rset_t **set)
309 {
310 if (*set != NULL) {
311 BFT_FREE((*set)->array);
312 BFT_FREE(*set);
313 }
314 }
315
316 /*----------------------------------------------------------------------------
317 * Check if we need to resize the current cs_join_rset_t structure and do
318 * it if necessary.
319 *
320 * parameters:
321 * set <-- pointer to pointer to the cs_join_rset_t structure to test
322 * test_size <-- target size
323 *---------------------------------------------------------------------------*/
324
325 void
cs_join_rset_resize(cs_join_rset_t ** set,cs_lnum_t test_size)326 cs_join_rset_resize(cs_join_rset_t **set,
327 cs_lnum_t test_size)
328 {
329 if (*set != NULL) {
330
331 if (test_size > 0) {
332
333 cs_join_rset_t *_set = *set;
334
335 if (test_size < _set->n_max_elts)
336 return;
337
338 if (_set->n_max_elts == 0)
339 _set->n_max_elts = test_size;
340 else if (test_size >= _set->n_max_elts) {
341 while (test_size >= _set->n_max_elts)
342 _set->n_max_elts *= 2; /* Double the list size */
343 }
344
345 BFT_REALLOC(_set->array, _set->n_max_elts, cs_lnum_t);
346 assert(test_size <= _set->n_max_elts);
347
348 }
349
350 }
351 else
352 *set = cs_join_rset_create(test_size);
353 }
354
355 /*----------------------------------------------------------------------------
356 * Create a new cs_join_eset_t structure.
357 *
358 * parameters:
359 * init_size <-- number of initial equivalences to allocate
360 *
361 * returns:
362 * a pointer to a new cs_join_eset_t structure
363 *---------------------------------------------------------------------------*/
364
365 cs_join_eset_t *
cs_join_eset_create(cs_lnum_t init_size)366 cs_join_eset_create(cs_lnum_t init_size)
367 {
368 cs_join_eset_t *new_set = NULL;
369
370 BFT_MALLOC(new_set, 1, cs_join_eset_t);
371
372 new_set->n_max_equiv = init_size; /* default value */
373 new_set->n_equiv = 0;
374
375 BFT_MALLOC(new_set->equiv_couple, 2*new_set->n_max_equiv, cs_lnum_t);
376
377 return new_set;
378 }
379
380 /*----------------------------------------------------------------------------
381 * Check if the requested size if allocated in the structure.
382 *
383 * Reallocate cs_join_eset_t structure if necessary.
384 *
385 * parameters:
386 * request_size <-- necessary size
387 * equiv_set <-> pointer to pointer to the cs_join_eset_t struct.
388 *---------------------------------------------------------------------------*/
389
390 void
cs_join_eset_check_size(cs_lnum_t request_size,cs_join_eset_t ** equiv_set)391 cs_join_eset_check_size(cs_lnum_t request_size,
392 cs_join_eset_t **equiv_set)
393 {
394 assert(equiv_set != NULL);
395
396 cs_join_eset_t *eset = *equiv_set;
397
398 if (eset == NULL)
399 eset = cs_join_eset_create(request_size);
400
401 if (request_size + 1 > eset->n_max_equiv) {
402
403 if (eset->n_max_equiv == 0)
404 eset->n_max_equiv = 2;
405
406 eset->n_max_equiv *= 2;
407
408 BFT_REALLOC(eset->equiv_couple, 2*eset->n_max_equiv, cs_lnum_t);
409
410 }
411
412 /* Ensure return value is set */
413
414 *equiv_set = eset;
415 }
416
417 /*----------------------------------------------------------------------------
418 * Destroy a cs_join_eset_t structure.
419 *
420 * parameters:
421 * equiv_set <-- pointer to pointer to the structure to destroy
422 *---------------------------------------------------------------------------*/
423
424 void
cs_join_eset_destroy(cs_join_eset_t ** equiv_set)425 cs_join_eset_destroy(cs_join_eset_t **equiv_set)
426 {
427 if (*equiv_set != NULL) {
428 BFT_FREE((*equiv_set)->equiv_couple);
429 BFT_FREE(*equiv_set);
430 }
431 }
432
433 /*----------------------------------------------------------------------------
434 * Clean a cs_join_eset_t structure.
435 *
436 * If necessary, create a new cs_join_eset_t structure with no redundancy.
437 *
438 * parameters:
439 * eset <-- pointer to pointer to the cs_join_eset_t structure to clean
440 *---------------------------------------------------------------------------*/
441
442 void
cs_join_eset_clean(cs_join_eset_t ** eset)443 cs_join_eset_clean(cs_join_eset_t **eset)
444 {
445 assert(eset != NULL);
446
447 int i;
448 cs_lnum_t prev, current;
449
450 cs_lnum_t count = 0;
451 cs_lnum_t *order = NULL;
452 cs_join_eset_t *new_eset = NULL;
453 cs_join_eset_t *_eset = *eset;
454
455 if (_eset == NULL)
456 return;
457
458 if (_eset->n_equiv == 1)
459 return;
460
461 BFT_MALLOC(order, _eset->n_equiv, cs_lnum_t);
462
463 if (_order_local_test_s2(_eset->equiv_couple,
464 _eset->n_equiv) == false) {
465
466 /* Order equiv_lst */
467
468 _order_local_s2(_eset->equiv_couple,
469 order,
470 _eset->n_equiv);
471
472 }
473 else
474 for (i = 0; i < _eset->n_equiv; i++)
475 order[i] = i;
476
477 /* Count the number of redundancies */
478
479 count = 0;
480
481 for (i = 1; i < _eset->n_equiv; i++) {
482
483 prev = order[i-1];
484 current = order[i];
485
486 if (_eset->equiv_couple[2*prev] == _eset->equiv_couple[2*current])
487 if (_eset->equiv_couple[2*prev+1] == _eset->equiv_couple[2*current+1])
488 count++;
489
490 } /* End of loop on equivalences */
491
492 new_eset = cs_join_eset_create(_eset->n_equiv - count);
493
494 new_eset->n_equiv = _eset->n_equiv - count;
495
496 if (new_eset->n_equiv > new_eset->n_max_equiv) {
497 new_eset->n_max_equiv = new_eset->n_equiv;
498 BFT_REALLOC(new_eset->equiv_couple, 2*new_eset->n_max_equiv, cs_lnum_t);
499 }
500
501 if (new_eset->n_equiv > 0) {
502
503 new_eset->equiv_couple[0] = _eset->equiv_couple[2*order[0]];
504 new_eset->equiv_couple[1] = _eset->equiv_couple[2*order[0]+1];
505 count = 1;
506
507 for (i = 1; i < _eset->n_equiv; i++) {
508
509 prev = order[i-1];
510 current = order[i];
511
512 if (_eset->equiv_couple[2*prev] != _eset->equiv_couple[2*current]) {
513 new_eset->equiv_couple[2*count] = _eset->equiv_couple[2*current];
514 new_eset->equiv_couple[2*count+1] = _eset->equiv_couple[2*current+1];
515 count++;
516 }
517 else {
518
519 if (_eset->equiv_couple[2*prev+1] != _eset->equiv_couple[2*current+1]) {
520 new_eset->equiv_couple[2*count] = _eset->equiv_couple[2*current];
521 new_eset->equiv_couple[2*count+1] = _eset->equiv_couple[2*current+1];
522 count++;
523 }
524
525 }
526
527 } /* End of loop on equivalences */
528
529 assert(count == new_eset->n_equiv);
530
531 }
532
533 *eset = new_eset;
534
535 /* Free memory */
536
537 cs_join_eset_destroy(&_eset);
538
539 BFT_FREE(order);
540 }
541
542 /*----------------------------------------------------------------------------
543 * Create a cs_join_gset_t structure (indexed list on global numbering)
544 *
545 * parameters:
546 * n_elts <-- number of elements composing the list
547 *
548 * returns:
549 * a new allocated pointer to a cs_join_gset_t structure.
550 *---------------------------------------------------------------------------*/
551
552 cs_join_gset_t *
cs_join_gset_create(cs_lnum_t n_elts)553 cs_join_gset_create(cs_lnum_t n_elts)
554 {
555 cs_lnum_t i;
556
557 cs_join_gset_t *new_set = NULL;
558
559 BFT_MALLOC(new_set, 1, cs_join_gset_t);
560
561 new_set->n_elts = n_elts;
562 new_set->n_g_elts = 0;
563
564 BFT_MALLOC(new_set->g_elts, n_elts, cs_gnum_t);
565 for (i = 0; i < n_elts; i++)
566 new_set->g_elts[i] = 0;
567
568 BFT_MALLOC(new_set->index, n_elts + 1, cs_lnum_t);
569 for (i = 0; i < n_elts + 1; i++)
570 new_set->index[i] = 0;
571
572 new_set->g_list = NULL;
573
574 return new_set;
575 }
576
577 /*----------------------------------------------------------------------------
578 * Build a cs_join_gset_t structure to store all the potential groups
579 * between elements.
580 *
581 * Values in g_elts are the tag values and values in g_list
582 * are position in tag array.
583 *
584 * parameters:
585 * n_elts <-- number of elements in tag array
586 * tag <-- tag array used to define a new cs_join_gset_t
587 *
588 * returns:
589 * a new allocated cs_join_gset_t structure
590 *---------------------------------------------------------------------------*/
591
592 cs_join_gset_t *
cs_join_gset_create_from_tag(cs_lnum_t n_elts,const cs_gnum_t tag[])593 cs_join_gset_create_from_tag(cs_lnum_t n_elts,
594 const cs_gnum_t tag[])
595 {
596 cs_lnum_t i, n_list_elts;
597 cs_gnum_t prev;
598
599 cs_lnum_t *order = NULL;
600 cs_join_gset_t *set = NULL;
601
602 if (n_elts == 0) {
603 set = cs_join_gset_create(n_elts);
604 return set;
605 }
606
607 /* Order tag */
608
609 assert(tag != NULL);
610
611 BFT_MALLOC(order, n_elts, cs_lnum_t);
612
613 cs_order_gnum_allocated(NULL, tag, order, n_elts);
614
615 /* Create a cs_join_gset_t structure to store the initial position of equiv.
616 element in tag */
617
618 prev = tag[order[0]];
619 n_list_elts = 1;
620
621 /* Count the number of elements which will compose the set->g_elts */
622
623 for (i = 1; i < n_elts; i++) {
624
625 cs_gnum_t cur = tag[order[i]];
626
627 if (prev != cur) {
628 n_list_elts++;
629 prev = cur;
630 }
631
632 }
633
634 set = cs_join_gset_create(n_list_elts);
635
636 if (n_list_elts > 0) {
637
638 cs_lnum_t shift;
639 cs_lnum_t count = 0;
640
641 /* Define the list of elements in set->g_elts and count the number of
642 associated entities */
643
644 prev = tag[order[0]];
645 set->g_elts[0] = prev;
646 set->index[1] += 1;
647 n_list_elts = 1;
648
649 for (i = 1; i < n_elts; i++) {
650
651 cs_gnum_t cur = tag[order[i]];
652
653 if (prev != cur) {
654 prev = cur;
655 set->g_elts[n_list_elts] = cur;
656 n_list_elts++;
657 set->index[n_list_elts] += 1;
658 }
659 else
660 set->index[n_list_elts] += 1;
661
662 }
663
664 /* Build index for the set */
665
666 for (i = 0; i < set->n_elts; i++)
667 set->index[i+1] += set->index[i];
668
669 /* Fill list */
670
671 BFT_MALLOC(set->g_list, set->index[set->n_elts], cs_gnum_t);
672
673 n_list_elts = 0;
674 prev = tag[order[0]];
675 set->g_list[0] = order[0];
676
677 for (i = 1; i < n_elts; i++) {
678
679 cs_lnum_t o_id = order[i];
680 cs_gnum_t cur = tag[o_id];
681
682 if (prev != cur) {
683 prev = cur;
684 count = 0;
685 n_list_elts++;
686 shift = set->index[n_list_elts];
687 set->g_list[shift] = o_id;
688 }
689 else {
690 count++;
691 shift = count + set->index[n_list_elts];
692 set->g_list[shift] = o_id;
693 }
694
695 }
696
697 } /* End if n_elts > 0 */
698
699 /* Free memory */
700
701 BFT_FREE(order);
702
703 /* Returns pointers */
704
705 return set;
706 }
707
708 /*----------------------------------------------------------------------------
709 * Create a new cs_join_gset_t which holds equivalences between elements of
710 * g_list in cs_join_gset_t.
711 *
712 * For a subset of equivalences, we store their initial value in the return
713 * cs_join_gset_t structure. A subset is defined if at least two elements
714 * are equivalent.
715 *
716 * The behavior of this function is near from cs_join_gset_create_from_tag
717 * but we don't store the position in init_array but its value in init_array.
718 *
719 * parameters:
720 * set <-- pointer to a cs_join_gset_t structure
721 * init_array <-- initial values of set->g_list
722 *
723 * returns:
724 * a new allocated cs_join_gset_t structure, or NULL if it would be empty
725 *---------------------------------------------------------------------------*/
726
727 cs_join_gset_t *
cs_join_gset_create_by_equiv(const cs_join_gset_t * set,const cs_gnum_t init_array[])728 cs_join_gset_create_by_equiv(const cs_join_gset_t *set,
729 const cs_gnum_t init_array[])
730 {
731 cs_lnum_t i, list_size, n_equiv_grp, count, shift, o_id;
732 cs_gnum_t prev, cur;
733
734 cs_lnum_t save_i = -1;
735 cs_lnum_t *order = NULL;
736 cs_gnum_t *couple_list = NULL;
737 cs_join_gset_t *equiv = NULL;
738
739 if (init_array == NULL)
740 return NULL;
741
742 assert(set != NULL);
743
744 list_size = set->index[set->n_elts];
745
746 /* Order set->g_list */
747
748 BFT_MALLOC(order, list_size, cs_lnum_t);
749 BFT_MALLOC(couple_list, 2*list_size, cs_gnum_t);
750
751 for (i = 0; i < list_size; i++) {
752 couple_list[2*i] = set->g_list[i];
753 couple_list[2*i+1] = init_array[i];
754 }
755
756 cs_order_gnum_allocated_s(NULL, couple_list, 2, order, list_size);
757
758 /* Create a cs_join_gset_t structure to store the initial value of equiv.
759 element in set->g_list */
760
761 /* Count the number of elements which will compose the equiv->g_elts */
762
763 n_equiv_grp = 0;
764 prev = set->g_list[order[0]];
765 count = 0;
766
767 for (i = 1; i < list_size; i++) {
768
769 cur = set->g_list[order[i]];
770
771 if (prev != cur) {
772 count = 0;
773 prev = cur;
774 }
775 else {
776 count++;
777 if (count == 1)
778 n_equiv_grp++;
779 }
780
781 }
782
783 equiv = cs_join_gset_create(n_equiv_grp);
784
785 if (n_equiv_grp > 0) {
786
787 /* Define the list of elements in equiv->g_list and count the number of
788 associated elements */
789
790 n_equiv_grp = 0;
791 prev = set->g_list[order[0]];
792 count = 0;
793
794 for (i = 1; i < list_size; i++) {
795
796 cur = set->g_list[order[i]];
797
798 if (prev != cur) {
799 count = 0;
800 prev = cur;
801 }
802 else {
803 count++;
804 if (count == 1) { /* Add this group */
805 equiv->g_elts[n_equiv_grp] = cur;
806 n_equiv_grp++;
807 equiv->index[n_equiv_grp] = 1;
808 }
809 else {
810 assert(count > 1);
811 equiv->index[n_equiv_grp] += 1;
812 }
813
814 } /* cur == prev */
815
816 } /* End of loop on list_size */
817
818 /* Build index for the set */
819
820 for (i = 0; i < equiv->n_elts; i++)
821 equiv->index[i+1] += equiv->index[i];
822
823 /* Fill list */
824
825 BFT_MALLOC(equiv->g_list, equiv->index[equiv->n_elts], cs_gnum_t);
826
827 n_equiv_grp = 0;
828 prev = set->g_list[order[0]] + 1;
829
830 for (i = 0; i < list_size; i++) {
831
832 o_id = order[i];
833 cur = set->g_list[o_id];
834
835 if (prev != cur) {
836 count = 0;
837 prev = cur;
838 save_i = o_id;
839 }
840 else {
841
842 if (count == 0)
843 n_equiv_grp++;
844
845 shift = count + equiv->index[n_equiv_grp-1];
846 if (cur != init_array[o_id])
847 equiv->g_list[shift] = init_array[o_id];
848 else
849 equiv->g_list[shift] = init_array[save_i];
850 count++;
851
852 } /* cur == prev */
853
854 } /* End of loop on list_size */
855
856 } /* End if n_elts > 0 */
857
858 /* Free memory */
859
860 BFT_FREE(couple_list);
861 BFT_FREE(order);
862
863 /* Returns pointer */
864
865 return equiv;
866 }
867
868 /*----------------------------------------------------------------------------
869 * Copy a cs_join_gset_t structure.
870 *
871 * parameters:
872 * src <-- pointer to the cs_join_gset_t structure to copy
873 *
874 * returns:
875 * a new allocated cs_join_gset_t structure.
876 *---------------------------------------------------------------------------*/
877
878 cs_join_gset_t *
cs_join_gset_copy(const cs_join_gset_t * src)879 cs_join_gset_copy(const cs_join_gset_t *src)
880 {
881 cs_lnum_t i;
882
883 cs_join_gset_t *copy = NULL;
884
885 if (src == NULL)
886 return copy;
887
888 copy = cs_join_gset_create(src->n_elts);
889
890 for (i = 0; i < src->n_elts; i++)
891 copy->g_elts[i] = src->g_elts[i];
892
893 for (i = 0; i < src->n_elts + 1; i++)
894 copy->index[i] = src->index[i];
895
896 BFT_MALLOC(copy->g_list, copy->index[copy->n_elts], cs_gnum_t);
897
898 for (i = 0; i < src->index[src->n_elts]; i++)
899 copy->g_list[i] = src->g_list[i];
900
901 return copy;
902 }
903
904 /*----------------------------------------------------------------------------
905 * Destroy a cs_join_gset_t structure.
906 *
907 * parameters:
908 * set <-- pointer to pointer to the cs_join_gset_t structure to destroy
909 *---------------------------------------------------------------------------*/
910
911 void
cs_join_gset_destroy(cs_join_gset_t ** set)912 cs_join_gset_destroy(cs_join_gset_t **set)
913 {
914 if (*set != NULL) {
915 BFT_FREE((*set)->index);
916 BFT_FREE((*set)->g_elts);
917 BFT_FREE((*set)->g_list);
918 BFT_FREE(*set);
919 }
920 }
921
922 /*----------------------------------------------------------------------------
923 * Sort a cs_join_gset_t structure according to the global numbering of
924 * the g_elts in cs_join_gset_t structure.
925 *
926 * parameters:
927 * set <-> pointer to the structure to order
928 *---------------------------------------------------------------------------*/
929
930 void
cs_join_gset_sort_elts(cs_join_gset_t * set)931 cs_join_gset_sort_elts(cs_join_gset_t *set)
932 {
933 int i, j, k, o_id, shift;
934 cs_lnum_t n_elts;
935
936 cs_lnum_t *new_index = NULL;
937 cs_lnum_t *order = NULL;
938 cs_gnum_t *tmp = NULL, *g_elts = NULL, *g_list = NULL;
939
940 if (set == NULL)
941 return;
942
943 g_elts = set->g_elts;
944 g_list = set->g_list;
945 n_elts = set->n_elts;
946
947 BFT_MALLOC(order, n_elts, cs_lnum_t);
948 BFT_MALLOC(tmp, n_elts, cs_gnum_t);
949 BFT_MALLOC(new_index, n_elts + 1, cs_lnum_t);
950
951 for (i = 0; i < n_elts; i++)
952 tmp[i] = g_elts[i];
953
954 /* Sort g_elts */
955
956 cs_order_gnum_allocated(NULL, g_elts, order, n_elts);
957
958 /* Reshape cs_join_gset_t according to the new ordering */
959
960 new_index[0] = 0;
961
962 for (i = 0; i < n_elts; i++) {
963
964 o_id = order[i];
965 g_elts[i] = tmp[o_id];
966 new_index[i+1] = new_index[i] + set->index[o_id+1] - set->index[o_id];
967
968 } /* End of loop on elements */
969
970 assert(new_index[n_elts] == set->index[n_elts]);
971
972 /* Define new g_list */
973
974 BFT_REALLOC(tmp, set->index[n_elts], cs_gnum_t);
975
976 for (i = 0; i < set->index[n_elts]; i++)
977 tmp[i] = g_list[i];
978
979 for (i = 0; i < n_elts; i++) {
980
981 o_id = order[i];
982 shift = new_index[i];
983
984 for (k = 0, j = set->index[o_id]; j < set->index[o_id+1]; j++, k++)
985 g_list[shift + k] = tmp[j];
986
987 } /* End of loop on elements */
988
989 /* Free memory */
990
991 BFT_FREE(set->index);
992 BFT_FREE(order);
993 BFT_FREE(tmp);
994
995 /* Return structure */
996
997 set->index = new_index;
998 set->g_elts = g_elts;
999 set->g_list = g_list;
1000 }
1001
1002 /*----------------------------------------------------------------------------
1003 * Sort each sub-list of the g_list array in a cs_join_gset_t structure.
1004 *
1005 * parameters:
1006 * p_set <-> pointer to the structure to sort
1007 *---------------------------------------------------------------------------*/
1008
1009 void
cs_join_gset_sort_sublist(cs_join_gset_t * set)1010 cs_join_gset_sort_sublist(cs_join_gset_t *set)
1011 {
1012 int i;
1013
1014 if (set == NULL)
1015 return;
1016
1017 /* Sort each sub-list */
1018
1019 for (i = 0; i < set->n_elts; i++)
1020 cs_sort_gnum_shell(set->index[i], set->index[i+1], set->g_list);
1021 }
1022
1023 /*----------------------------------------------------------------------------
1024 * Invert a cs_join_gset_t structure.
1025 *
1026 * parameters:
1027 * set <-- pointer to the cs_join_gset_t structure to work with
1028 *
1029 * returns:
1030 * the new allocated and inverted set structure
1031 *---------------------------------------------------------------------------*/
1032
1033 cs_join_gset_t *
cs_join_gset_invert(const cs_join_gset_t * set)1034 cs_join_gset_invert(const cs_join_gset_t *set)
1035 {
1036 int i, j, o_id, shift, elt_id;
1037 cs_gnum_t prev, cur;
1038
1039 cs_lnum_t list_size = 0, n_elts = 0;
1040 cs_lnum_t *count = NULL;
1041 cs_lnum_t *order = NULL;
1042 cs_join_gset_t *invert_set = NULL;
1043
1044 if (set == NULL)
1045 return invert_set;
1046
1047 /* Order g_list to count the number of different entities */
1048
1049 list_size = set->index[set->n_elts];
1050
1051 if (list_size == 0)
1052 return cs_join_gset_create(list_size);
1053
1054 BFT_MALLOC(order, list_size, cs_lnum_t);
1055
1056 cs_order_gnum_allocated(NULL, set->g_list, order, list_size);
1057
1058 /* Count the number of elements */
1059
1060 prev = set->g_list[order[0]] + 1;
1061
1062 for (i = 0; i < list_size; i++) {
1063
1064 o_id = order[i];
1065 cur = set->g_list[o_id];
1066
1067 if (prev != cur) {
1068 prev = cur;
1069 n_elts++;
1070 }
1071
1072 }
1073
1074 invert_set = cs_join_gset_create(n_elts);
1075
1076 /* Fill g_elts for the inverted set */
1077
1078 n_elts = 0;
1079 prev = set->g_list[order[0]] + 1;
1080
1081 for (i = 0; i < list_size; i++) {
1082
1083 o_id = order[i];
1084 cur = set->g_list[o_id];
1085
1086 if (prev != cur) {
1087 prev = cur;
1088 invert_set->g_elts[n_elts] = cur;
1089 n_elts++;
1090 }
1091
1092 }
1093
1094 BFT_FREE(order);
1095
1096 /* Define an index for the inverted set */
1097
1098 for (i = 0; i < set->n_elts; i++) {
1099 for (j = set->index[i]; j < set->index[i+1]; j++) {
1100
1101 elt_id = cs_search_g_binary(invert_set->n_elts,
1102 set->g_list[j],
1103 invert_set->g_elts);
1104
1105 if (elt_id == -1)
1106 bft_error(__FILE__, __LINE__, 0,
1107 _(" Fail to build an inverted cs_join_gset_t structure.\n"
1108 " Cannot find %llu in element list.\n"),
1109 (unsigned long long)(set->g_list[j]));
1110
1111 invert_set->index[elt_id+1] += 1;
1112
1113 }
1114 } /* End of loop on elements of list */
1115
1116 for (i = 0; i < invert_set->n_elts; i++)
1117 invert_set->index[i+1] += invert_set->index[i];
1118
1119 BFT_MALLOC(invert_set->g_list,
1120 invert_set->index[invert_set->n_elts],
1121 cs_gnum_t);
1122
1123 /* Define invert_set->g_list */
1124
1125 BFT_MALLOC(count, invert_set->n_elts, cs_lnum_t);
1126
1127 for (i = 0; i < invert_set->n_elts; i++)
1128 count[i] = 0;
1129
1130 for (i = 0; i < set->n_elts; i++) {
1131 for (j = set->index[i]; j < set->index[i+1]; j++) {
1132
1133 elt_id = cs_search_g_binary(invert_set->n_elts,
1134 set->g_list[j],
1135 invert_set->g_elts);
1136
1137 shift = count[elt_id] + invert_set->index[elt_id];
1138 invert_set->g_list[shift] = set->g_elts[i];
1139 count[elt_id] += 1;
1140
1141 }
1142
1143 } /* End of loop on elements of list */
1144
1145 BFT_FREE(count);
1146
1147 return invert_set;
1148 }
1149
1150 /*----------------------------------------------------------------------------
1151 * Delete redudancies in a cs_join_gset_t structure.
1152 *
1153 * Output set has an ordered sub-list for each element in set.
1154 *
1155 * parameters:
1156 * set <-> pointer to the structure to clean
1157 *---------------------------------------------------------------------------*/
1158
1159 void
cs_join_gset_clean(cs_join_gset_t * set)1160 cs_join_gset_clean(cs_join_gset_t *set)
1161 {
1162 int i, j, l, r, save, n_elts;
1163
1164 int shift = 0;
1165 cs_gnum_t *g_list = NULL;
1166
1167 if (set == NULL)
1168 return;
1169
1170 g_list = set->g_list;
1171 n_elts = set->n_elts;
1172
1173 /* Sort g_list for each element in index */
1174
1175 cs_join_gset_sort_sublist(set);
1176
1177 /* Define a new index without redundant elements */
1178
1179 save = set->index[0];
1180
1181 for (i = 0; i < n_elts; i++) {
1182
1183 l = save;
1184 r = set->index[i+1];
1185
1186 if (r - l > 0) {
1187
1188 g_list[shift++] = g_list[l];
1189
1190 for (j = l + 1; j < r; j++)
1191 if (g_list[j] != g_list[j-1])
1192 g_list[shift++] = g_list[j];
1193
1194 }
1195
1196 save = r;
1197 set->index[i+1] = shift;
1198
1199 } /* End of loop on elements */
1200 }
1201
1202 /*----------------------------------------------------------------------------
1203 * Delete redudancies in g_list array of a cs_join_gset_t structure.
1204 *
1205 * parameters:
1206 * set <-> pointer to the structure to clean
1207 * linked_array <-> array for which redundancies are scanned
1208 *---------------------------------------------------------------------------*/
1209
1210 void
cs_join_gset_clean_from_array(cs_join_gset_t * set,cs_gnum_t linked_array[])1211 cs_join_gset_clean_from_array(cs_join_gset_t *set,
1212 cs_gnum_t linked_array[])
1213 {
1214 int i, j, l, r;
1215 cs_lnum_t n_elts;
1216
1217 int shift = 0;
1218 cs_lnum_t *new_index = NULL;
1219 cs_gnum_t *g_list = NULL;
1220
1221 if (set == NULL)
1222 return;
1223
1224 if (linked_array == NULL)
1225 return;
1226
1227 g_list = set->g_list;
1228 n_elts = set->n_elts;
1229
1230 /* Sort linked_array and apply change to g_list for each element in index */
1231
1232 for (i = 0; i < n_elts; i++)
1233 _coupled_adapted_gnum_shellsort(set->index[i],
1234 set->index[i+1],
1235 linked_array,
1236 g_list);
1237
1238 /* Define a new index without redundant elements */
1239
1240 BFT_MALLOC(new_index, n_elts + 1, cs_lnum_t);
1241 new_index[0] = 0;
1242
1243 for (i = 0; i < n_elts; i++) {
1244
1245 l = set->index[i];
1246 r = set->index[i+1];
1247
1248 if (r - l > 0) {
1249
1250 g_list[shift] = g_list[l];
1251 shift++;
1252
1253 for (j = l + 1; j < r; j++) {
1254
1255 if (linked_array[j] != linked_array[j-1]) {
1256 g_list[shift] = g_list[j];
1257 shift++;
1258 }
1259 assert(g_list[shift-1] <= g_list[j]);
1260
1261 }
1262 new_index[i+1] = shift;
1263
1264 }
1265 else { /* No match for this element */
1266
1267 new_index[i+1] = new_index[i];
1268
1269 }
1270
1271 } /* End of loop on elements */
1272
1273 /* Reshape cs_join_gset_t structure */
1274
1275 BFT_REALLOC(g_list, new_index[n_elts], cs_gnum_t);
1276 BFT_FREE(set->index);
1277
1278 set->index = new_index;
1279 set->g_list = g_list;
1280 }
1281
1282 /*----------------------------------------------------------------------------
1283 * Concatenate the two g_elts and g_list arrays.
1284 *
1285 * Order the new concatenated array and delete redundant elements.
1286 * We get a single ordered array.
1287 *
1288 * parameters:
1289 * set <-- pointer to the structure to work with
1290 * n_elts --> number of elements in the new set
1291 * new_array --> pointer to the new created array
1292 *---------------------------------------------------------------------------*/
1293
1294 void
cs_join_gset_single_order(const cs_join_gset_t * set,cs_lnum_t * n_elts,cs_gnum_t * new_array[])1295 cs_join_gset_single_order(const cs_join_gset_t *set,
1296 cs_lnum_t *n_elts,
1297 cs_gnum_t *new_array[])
1298 {
1299 cs_lnum_t _n_elts = 0;
1300 cs_gnum_t *_new_array = NULL;
1301
1302 *n_elts = _n_elts;
1303 *new_array = _new_array;
1304
1305 if (set == NULL) /* Nothing to do */
1306 return;
1307
1308 _n_elts = set->n_elts;
1309
1310 if (_n_elts > 0) {
1311
1312 cs_lnum_t i, shift;
1313 cs_gnum_t prev;
1314
1315 cs_lnum_t *order = NULL;
1316 cs_gnum_t *elt_list = NULL;
1317
1318 _n_elts += set->index[_n_elts]; /* Add number of elements in g_list */
1319
1320 BFT_MALLOC(elt_list, _n_elts, cs_gnum_t);
1321
1322 for (i = 0; i < set->n_elts; i++)
1323 elt_list[i] = set->g_elts[i];
1324
1325 shift = set->n_elts;
1326 for (i = 0; i < set->index[set->n_elts]; i++)
1327 elt_list[shift + i] = set->g_list[i];
1328
1329 /* Define an ordered list of elements */
1330
1331 BFT_MALLOC(_new_array, _n_elts, cs_gnum_t);
1332 BFT_MALLOC(order, _n_elts, cs_lnum_t);
1333
1334 cs_order_gnum_allocated(NULL, elt_list, order, _n_elts);
1335
1336 for (i = 0; i < _n_elts; i++)
1337 _new_array[i] = elt_list[order[i]];
1338
1339 /* Delete redundant elements */
1340
1341 shift = 0;
1342 prev = _new_array[0] + 1;
1343
1344 for (i = 0; i < _n_elts; i++) {
1345
1346 if (prev != _new_array[i]) {
1347
1348 _new_array[shift] = _new_array[i];
1349 prev = _new_array[i];
1350 shift++;
1351
1352 }
1353
1354 }
1355 _n_elts = shift; /* Real number of elements without redundancy */
1356
1357 /* Memory management */
1358
1359 BFT_FREE(order);
1360 BFT_FREE(elt_list);
1361 BFT_REALLOC(_new_array, _n_elts, cs_gnum_t);
1362
1363 } /* End if n_elts > 0 */
1364
1365 /* Set return pointers */
1366
1367 *n_elts = _n_elts;
1368 *new_array = _new_array;
1369 }
1370
1371 /*----------------------------------------------------------------------------
1372 * Compress a g_list such as for each element "e" in g_elts:
1373 * - there is no redundancy for the linked elements of set->g_list
1374 * - there is no element in set->g_list < e except if this element is not
1375 * present in g_elts
1376 *
1377 * g_list and g_elts must be ordered before calling this function
1378 *
1379 * parameters:
1380 * set <-> pointer to the structure to work with
1381 *---------------------------------------------------------------------------*/
1382
1383 void
cs_join_gset_compress(cs_join_gset_t * set)1384 cs_join_gset_compress(cs_join_gset_t *set)
1385 {
1386 cs_lnum_t i, j, start, end, save, shift;
1387 cs_gnum_t cur;
1388
1389 if (set == NULL)
1390 return;
1391
1392 if (set->n_elts == 0)
1393 return;
1394
1395 shift = 0;
1396 save = set->index[0];
1397
1398 for (i = 0; i < set->n_elts; i++) {
1399
1400 cur = set->g_elts[i];
1401 start = save;
1402 end = set->index[i+1];
1403
1404 if (end - start > 0) {
1405
1406 /* Sub-lists must be ordered */
1407
1408 if (cur < set->g_list[start])
1409 set->g_list[shift++] = set->g_list[start];
1410 else if (cur > set->g_list[start]) {
1411
1412 int id = cs_search_g_binary(i+1,
1413 set->g_list[start],
1414 set->g_elts);
1415
1416 if (id == -1) /* Not found. Keep it. */
1417 set->g_list[shift++] = set->g_list[start];
1418
1419 }
1420
1421 for (j = start + 1; j < end; j++) {
1422
1423 if (cur < set->g_list[j]) {
1424 if (set->g_list[j-1] != set->g_list[j])
1425 set->g_list[shift++] = set->g_list[j];
1426 }
1427 else if (cur > set->g_list[j]) {
1428
1429 int id = cs_search_g_binary(i+1,
1430 set->g_list[j],
1431 set->g_elts);
1432
1433 if (id == -1) /* Not found. Keep it. */
1434 set->g_list[shift++] = set->g_list[j];
1435
1436 }
1437
1438 } /* End of loop on sub-elements */
1439
1440 } /* end - start > 0 */
1441
1442 save = end;
1443 set->index[i+1] = shift;
1444
1445 } /* End of loop on elements in g_elts */
1446
1447 /* Reshape cs_join_gset_t structure if necessary */
1448
1449 if (save != set->index[set->n_elts]) {
1450 assert(save > set->index[set->n_elts]);
1451 BFT_REALLOC(set->g_list, set->index[set->n_elts], cs_gnum_t);
1452 }
1453
1454 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1455 cs_join_gset_dump(NULL, set);
1456 #endif
1457 }
1458
1459 /*----------------------------------------------------------------------------
1460 * Delete redundancies in set->g_elts.
1461 *
1462 * Merge sub-arrays associated to a common set->g_elts[i].
1463 *
1464 * parameters:
1465 * set <-- pointer to the structure to work with
1466 * order_tag <-- 0: set->g_elts is not ordered, 1: ordered
1467 *---------------------------------------------------------------------------*/
1468
1469 void
cs_join_gset_merge_elts(cs_join_gset_t * set,int order_tag)1470 cs_join_gset_merge_elts(cs_join_gset_t *set,
1471 int order_tag)
1472 {
1473 cs_lnum_t i, save, start, end, n_init_elts, n_sub_elts;
1474 cs_gnum_t prev, cur;
1475
1476 if (set == NULL)
1477 return;
1478
1479 n_init_elts = set->n_elts;
1480
1481 if (n_init_elts < 2)
1482 return;
1483
1484 assert(order_tag == 0 || order_tag == 1);
1485
1486 /* Delete redundancies in g_elts. Merge sub-lists associated to common
1487 g_elts */
1488
1489 if (order_tag == 0)
1490 cs_join_gset_sort_elts(set);
1491
1492 set->n_elts = 0; /* Reset and will be redefinied */
1493 prev = set->g_elts[0] + 1; /* Force prev to be different from g_elts[0] */
1494 save = set->index[0];
1495
1496 for (i = 0; i < n_init_elts; i++) {
1497
1498 cur = set->g_elts[i];
1499 start = save;
1500 end = set->index[i+1];
1501 save = end;
1502 n_sub_elts = end - start;
1503
1504 if (prev != cur) {
1505
1506 set->g_elts[set->n_elts] = cur;
1507 set->n_elts += 1;
1508 set->index[set->n_elts] = n_sub_elts;
1509 prev = cur;
1510
1511 }
1512 else {
1513
1514 set->index[set->n_elts] += n_sub_elts;
1515
1516 } /* prev != next */
1517
1518 } /* Loop on elements of set->g_elts */
1519
1520 /* Get the new index */
1521
1522 for (i = 0; i < set->n_elts; i++)
1523 set->index[i+1] += set->index[i];
1524
1525 /* Reshape cs_join_gset_t structure if necessary */
1526
1527 if (n_init_elts != set->n_elts) {
1528
1529 assert(n_init_elts > set->n_elts);
1530
1531 BFT_REALLOC(set->g_elts, set->n_elts, cs_gnum_t);
1532 BFT_REALLOC(set->index, set->n_elts + 1, cs_lnum_t);
1533 BFT_REALLOC(set->g_list, set->index[set->n_elts], cs_gnum_t);
1534
1535 }
1536
1537 #if 0 && defined(DEBUG) && !defined(NDEBUG)
1538 cs_join_gset_dump(NULL, set);
1539 #endif
1540 }
1541
1542 #if defined(HAVE_MPI)
1543
1544 /*----------------------------------------------------------------------------
1545 * Synchronize a cs_join_gset_t structure and distribute the resulting set
1546 * over the rank by block
1547 *
1548 * parameters:
1549 * max_gnum <-- max global number in global element numbering
1550 * loc_set <-- pointer to the local structure to work with
1551 * comm <-- mpi_comm on which synchro. and distribution take place
1552 *
1553 * returns:
1554 * a synchronized and distributed cs_join_gset_t structure.
1555 *---------------------------------------------------------------------------*/
1556
1557 cs_join_gset_t *
cs_join_gset_block_sync(cs_gnum_t max_gnum,cs_join_gset_t * loc_set,MPI_Comm comm)1558 cs_join_gset_block_sync(cs_gnum_t max_gnum,
1559 cs_join_gset_t *loc_set,
1560 MPI_Comm comm)
1561 {
1562 cs_join_gset_t *sync_set = NULL;
1563
1564 int local_rank, n_ranks;
1565
1566 assert(loc_set != NULL);
1567
1568 if (max_gnum == 0)
1569 return sync_set;
1570
1571 MPI_Comm_rank(comm, &local_rank);
1572 MPI_Comm_size(comm, &n_ranks);
1573
1574 cs_block_dist_info_t bi
1575 = cs_block_dist_compute_sizes(local_rank,
1576 n_ranks,
1577 1,
1578 0,
1579 max_gnum);
1580
1581 cs_lnum_t block_size = 0;
1582 if (bi.gnum_range[1] > bi.gnum_range[0])
1583 block_size = bi.gnum_range[1] - bi.gnum_range[0];
1584
1585 cs_all_to_all_t *
1586 d = cs_all_to_all_create_from_block(loc_set->n_elts,
1587 0, /* flags */
1588 loc_set->g_elts,
1589 bi,
1590 comm);
1591
1592 /* Synchronize list definition for each global element */
1593
1594 cs_lnum_t *p_index;
1595 BFT_MALLOC(p_index, loc_set->n_elts+1, cs_lnum_t);
1596 cs_gnum_t *p_buffer;
1597 BFT_MALLOC(p_buffer,
1598 loc_set->index[loc_set->n_elts] + loc_set->n_elts,
1599 cs_gnum_t);
1600
1601 p_index[0] = 0;
1602 for (cs_lnum_t i = 0; i < loc_set->n_elts; i++) {
1603
1604 cs_lnum_t shift = p_index[i];
1605 p_buffer[shift++] = loc_set->g_elts[i];
1606 cs_lnum_t s_id = loc_set->index[i], e_id = loc_set->index[i+1];
1607 for (cs_lnum_t j = s_id; j < e_id; j++)
1608 p_buffer[shift++] = loc_set->g_list[j];
1609
1610 p_index[i+1] = shift;
1611
1612 }
1613
1614 cs_lnum_t *r_index
1615 = cs_all_to_all_copy_index(d,
1616 false, /* reverse */
1617 p_index,
1618 NULL);
1619
1620 cs_gnum_t *r_buffer
1621 = cs_all_to_all_copy_indexed(d,
1622 CS_GNUM_TYPE,
1623 false, /* reverse */
1624 p_index,
1625 p_buffer,
1626 r_index,
1627 NULL);
1628
1629 BFT_FREE(p_index);
1630 BFT_FREE(p_buffer);
1631
1632 cs_lnum_t n_r_elts = cs_all_to_all_n_elts_dest(d);
1633
1634 cs_all_to_all_destroy(&d);
1635
1636 /* Define sync_set: a distributed cs_join_gset_t structure which
1637 synchronizes data over the ranks */
1638
1639 sync_set = cs_join_gset_create(block_size);
1640
1641 for (cs_lnum_t i = 0; i < sync_set->n_elts; i++) {
1642 sync_set->g_elts[i] = bi.gnum_range[0] + (cs_gnum_t)i;
1643 }
1644
1645 /* Build index */
1646
1647 for (cs_lnum_t i = 0; i < n_r_elts; i++) {
1648 cs_lnum_t j = r_buffer[r_index[i]] - bi.gnum_range[0];
1649 sync_set->index[j+1] += r_index[i+1] - r_index[i] - 1;
1650 }
1651
1652 for (cs_lnum_t i = 0; i < sync_set->n_elts; i++) {
1653 sync_set->index[i+1] += sync_set->index[i];
1654 }
1655
1656 BFT_MALLOC(sync_set->g_list,
1657 sync_set->index[sync_set->n_elts],
1658 cs_gnum_t);
1659
1660 /* Now build set */
1661
1662 cs_lnum_t *count;
1663 BFT_MALLOC(count, sync_set->n_elts, cs_lnum_t);
1664 for (cs_lnum_t i = 0; i < sync_set->n_elts; i++)
1665 count[i] = 0;
1666
1667 for (cs_lnum_t i = 0; i < n_r_elts; i++) {
1668 cs_lnum_t r_shift = r_index[i];
1669 cs_lnum_t j = r_buffer[r_shift] - bi.gnum_range[0];
1670 cs_lnum_t w_shift = sync_set->index[j] + count[j];
1671
1672 cs_lnum_t n_sub_elts = r_index[i+1] - r_index[i] - 1;
1673 for (cs_lnum_t k = 0; k < n_sub_elts; k++)
1674 sync_set->g_list[w_shift + k] = r_buffer[r_shift + 1 + k];
1675 count[j] += n_sub_elts;
1676 }
1677
1678 BFT_FREE(count);
1679 BFT_FREE(r_buffer);
1680 BFT_FREE(r_index);
1681
1682 /* Return pointer */
1683
1684 cs_join_gset_clean(sync_set);
1685
1686 return sync_set;
1687 }
1688
1689 /*----------------------------------------------------------------------------
1690 * Update a local cs_join_gset_t structure from a distributed and
1691 * synchronized cs_join_gset_t structure.
1692 *
1693 * Loc_set should not have redundant elements.
1694 *
1695 * parameters:
1696 * max_gnum <-- max global number in global element numbering
1697 * sync_set <-- pointer to the structure which holds a synchronized block
1698 * loc_set <-> pointer to a local structure holding elements to update
1699 * comm <-- comm on which synchronization and distribution take place
1700 *---------------------------------------------------------------------------*/
1701
1702 void
cs_join_gset_block_update(cs_gnum_t max_gnum,const cs_join_gset_t * sync_set,cs_join_gset_t * loc_set,MPI_Comm comm)1703 cs_join_gset_block_update(cs_gnum_t max_gnum,
1704 const cs_join_gset_t *sync_set,
1705 cs_join_gset_t *loc_set,
1706 MPI_Comm comm)
1707 {
1708 int local_rank, n_ranks;
1709
1710 if (max_gnum == 0)
1711 return;
1712
1713 /* Sanity checks */
1714
1715 assert(loc_set != NULL);
1716 assert(sync_set != NULL);
1717
1718 /* Build a cs_join_block_info_t structure */
1719
1720 MPI_Comm_rank(comm, &local_rank);
1721 MPI_Comm_size(comm, &n_ranks);
1722
1723 cs_block_dist_info_t bi = cs_block_dist_compute_sizes(local_rank,
1724 n_ranks,
1725 1,
1726 0,
1727 max_gnum);
1728
1729 cs_all_to_all_t
1730 *d = cs_all_to_all_create_from_block(loc_set->n_elts,
1731 0, /* flags */
1732 loc_set->g_elts,
1733 bi,
1734 comm);
1735
1736 cs_gnum_t *wanted_elts = cs_all_to_all_copy_array(d,
1737 CS_GNUM_TYPE,
1738 1,
1739 false, /* reverse */
1740 loc_set->g_elts,
1741 NULL);
1742
1743 cs_lnum_t n_recv_elts = cs_all_to_all_n_elts_dest(d);
1744
1745 /* Get a synchronized list definition for each global element */
1746
1747 /* Send new list definition holding by sync_set to ranks which have
1748 request it */
1749
1750 cs_lnum_t *block_index;
1751 BFT_MALLOC(block_index, n_recv_elts+1, cs_lnum_t);
1752
1753 block_index[0] = 0;
1754 for (cs_lnum_t i = 0; i < n_recv_elts; i++) {
1755 cs_lnum_t block_id = wanted_elts[i] - bi.gnum_range[0];
1756 cs_lnum_t n_sub_elts = sync_set->index[block_id+1]
1757 - sync_set->index[block_id];
1758 block_index[i+1] = block_index[i] + n_sub_elts;
1759 }
1760
1761 cs_all_to_all_copy_index(d,
1762 true, /* reverse */
1763 block_index,
1764 loc_set->index);
1765
1766 cs_gnum_t *block_tuples;
1767 BFT_MALLOC(block_tuples, block_index[n_recv_elts], cs_gnum_t);
1768
1769 for (cs_lnum_t i = 0, shift = 0; i < n_recv_elts; i++) {
1770
1771 cs_lnum_t block_id = wanted_elts[i] - bi.gnum_range[0];
1772
1773 cs_lnum_t s_id = sync_set->index[block_id];
1774 cs_lnum_t e_id = sync_set->index[block_id+1];
1775 cs_lnum_t n_sub_elts = e_id - s_id;
1776
1777 for (cs_lnum_t j = s_id, k = 0; j < e_id; j++, k++)
1778 block_tuples[shift+k] = sync_set->g_list[j];
1779
1780 shift += n_sub_elts;
1781
1782 }
1783
1784 BFT_FREE(wanted_elts);
1785
1786 /* Re-initialize loc_set */
1787
1788 BFT_FREE(loc_set->g_list);
1789
1790 loc_set->g_list = cs_all_to_all_copy_indexed(d,
1791 CS_GNUM_TYPE,
1792 true, /* reverse */
1793 block_index,
1794 block_tuples,
1795 loc_set->index,
1796 NULL);
1797
1798 /* Free some memory */
1799
1800 cs_all_to_all_destroy(&d);
1801
1802 BFT_FREE(block_index);
1803 BFT_FREE(block_tuples);
1804 }
1805
1806 #endif /* HAVE_MPI */
1807
1808 /*----------------------------------------------------------------------------
1809 * Dump an array (int or double).
1810 *
1811 * This function is called according to the verbosity.
1812 *
1813 * parameters:
1814 * f <-- handle to output file
1815 * type <-- type of the array to display
1816 * header <-- header to display in front of the array
1817 * n_elts <-- number of elements to display
1818 * array <-- array to display
1819 *---------------------------------------------------------------------------*/
1820
1821 void
cs_join_dump_array(FILE * f,const char * type,const char * header,int n_elts,const void * array)1822 cs_join_dump_array(FILE *f,
1823 const char *type,
1824 const char *header,
1825 int n_elts,
1826 const void *array)
1827 {
1828 int i;
1829
1830 fprintf(f, " %s: ", header);
1831
1832 if (!strncmp(type, "int", strlen("int"))) { /* "int" array */
1833
1834 const int *i_array = array;
1835
1836 for (i = 0; i < n_elts; i++)
1837 fprintf(f, " %8d", i_array[i]);
1838
1839 }
1840 else if (!strncmp(type, "bool", strlen("bool"))) { /* "boolean" array */
1841
1842 const bool *b_array = array;
1843
1844 for (i = 0; i < n_elts; i++)
1845 if (b_array[i] == true)
1846 fprintf(f, " T");
1847 else {
1848 assert(b_array[i] == false);
1849 fprintf(f, " F");
1850 }
1851 }
1852 else if (!strncmp(type, "double", strlen("double"))) { /* "double" array */
1853
1854 const double *d_array = array;
1855
1856 for (i = 0; i < n_elts; i++)
1857 fprintf(f, " %10.8e", d_array[i]);
1858
1859 }
1860 else if (!strncmp(type, "gnum", strlen("gnum"))) { /* "gnum" array */
1861
1862 const cs_gnum_t *u_array = array;
1863
1864 for (i = 0; i < n_elts; i++)
1865 fprintf(f, " %9llu", (unsigned long long)u_array[i]);
1866
1867 }
1868 else
1869 bft_error(__FILE__, __LINE__, 0,
1870 " Unexpected type (%s) to display in the current dump.\n",
1871 type);
1872
1873 fprintf(f, "\n");
1874 }
1875
1876 /*----------------------------------------------------------------------------
1877 * Dump a cs_join_gset_t structure.
1878 *
1879 * parameters:
1880 * f <-- handle to output file
1881 * set <-- pointer to the cs_join_gset_t structure to dump
1882 *---------------------------------------------------------------------------*/
1883
1884 void
cs_join_gset_dump(FILE * f,const cs_join_gset_t * set)1885 cs_join_gset_dump(FILE *f,
1886 const cs_join_gset_t *set)
1887 {
1888 int i, j;
1889
1890 if (set == NULL)
1891 return;
1892
1893 if (f == NULL)
1894 f = stdout;
1895
1896 fprintf(f, "\nDump cs_join_gset_t structure: %p\n", (const void *)set);
1897 fprintf(f, "number of elements: %10ld\n", (long)set->n_elts);
1898 fprintf(f, "size of the list : %10ld\n\n", (long)set->index[set->n_elts]);
1899
1900 for (i = 0; i < set->n_elts; i++) {
1901
1902 int s = set->index[i], e = set->index[i+1];
1903 int n_matches = e-s;
1904 int n_loops = n_matches/10;
1905
1906 fprintf(f, "Global num: %8llu | subsize: %3d |",
1907 (unsigned long long)set->g_elts[i], n_matches);
1908
1909 for (j = 0; j < n_loops; j++) {
1910 if (j == 0)
1911 fprintf(f,
1912 "%8llu %8llu %8llu %8llu %8llu "
1913 "%8llu %8llu %8llu %8llu %8llu\n",
1914 (unsigned long long)set->g_list[s+ 10*j],
1915 (unsigned long long)set->g_list[s+ 10*j + 1],
1916 (unsigned long long)set->g_list[s+ 10*j + 2],
1917 (unsigned long long)set->g_list[s+ 10*j + 3],
1918 (unsigned long long)set->g_list[s+ 10*j + 4],
1919 (unsigned long long)set->g_list[s+ 10*j + 5],
1920 (unsigned long long)set->g_list[s+ 10*j + 6],
1921 (unsigned long long)set->g_list[s+ 10*j + 7],
1922 (unsigned long long)set->g_list[s+ 10*j + 8],
1923 (unsigned long long)set->g_list[s+ 10*j + 9]);
1924 else
1925 fprintf(f, " "
1926 "%8llu %8llu %8llu %8llu %8llu "
1927 "%8llu %8llu %8llu %8llu %8llu\n",
1928 (unsigned long long)set->g_list[s+ 10*j],
1929 (unsigned long long)set->g_list[s+ 10*j + 1],
1930 (unsigned long long)set->g_list[s+ 10*j + 2],
1931 (unsigned long long)set->g_list[s+ 10*j + 3],
1932 (unsigned long long)set->g_list[s+ 10*j + 4],
1933 (unsigned long long)set->g_list[s+ 10*j + 5],
1934 (unsigned long long)set->g_list[s+ 10*j + 6],
1935 (unsigned long long)set->g_list[s+ 10*j + 7],
1936 (unsigned long long)set->g_list[s+ 10*j + 8],
1937 (unsigned long long)set->g_list[s+ 10*j + 9]);
1938 }
1939
1940 if (e - s+10*n_loops > 0) {
1941 for (j = s + 10*n_loops; j < e; j++) {
1942 if (j == s + 10*n_loops && n_loops > 0)
1943 fprintf(f, " ");
1944 fprintf(f, "%8llu ", (unsigned long long)set->g_list[j]);
1945 }
1946 fprintf(f, "\n");
1947 }
1948
1949 if (n_matches == 0)
1950 fprintf(f, "\n");
1951
1952 } /* End of loop on boxes */
1953
1954 fflush(f);
1955 }
1956
1957 /*---------------------------------------------------------------------------*/
1958
1959 END_C_DECLS
1960