BELFEM 0.9.3
Berkeley Lab Finite Element Framework
Loading...
Searching...
No Matches
commtools.hpp
Go to the documentation of this file.
1/*
2 * BELFEM -- The Berkeley Lab Finite Element Framework
3 * Copyright (c) 2026, The Regents of the University of California,
4 * through Lawrence Berkeley National Laboratory (subject to receipt of any required
5 * approvals from the U.S. Dept. of Energy). All rights reserved.
6 *
7 * Developers: Christian Messe, Gregory Giard
8 *
9 * See the top-level LICENSE file for the complete license and disclaimer.
10 */
11
12#ifndef COMMTOOLS_HPP
13#define COMMTOOLS_HPP
14
15#include <limits>
16#include <type_traits>
17#include <string>
18
19#include "typedefs.hpp"
20#include "assert.hpp"
21#include "cl_Vector.hpp"
22#include "cl_Matrix.hpp"
23#include "cl_Communicator.hpp"
24#include "commtypes.hpp"
25
26extern belfem::Communicator gComm;
27
28namespace belfem
29{
30
31 constexpr int gMaxCommChunkLength = 64 * 1024;
32
33//==============================================================================
34// HELPERS TO AVOID AMBIGUITIES
35//==============================================================================
36
37 template<typename T>
38 struct is_scalar : std::true_type {};
39
40 template<typename T>
41 struct is_scalar<Vector<T>> : std::false_type {};
42
43 template<typename T>
44 struct is_scalar<Matrix<T>> : std::false_type {};
45
46//==============================================================================
47// UTILITIES
48//==============================================================================
49
50 proc_t
51 comm_size();
52
53//------------------------------------------------------------------------------
54
55 proc_t
56 comm_rank();
57
58//------------------------------------------------------------------------------
59
62 void
63 comm_check( const int aErrorCode );
64
65//------------------------------------------------------------------------------
66
68 void
69 comm_barrier();
70
71//------------------------------------------------------------------------------
72
100 void
101 comm_drain_check( const char * aLabel );
102
103//------------------------------------------------------------------------------
104
109 int
110 comm_tag( const proc_t aSource, const proc_t aTarget );
111
112//------------------------------------------------------------------------------
113
117 Cell< int >
118 comm_split( const index_t aLength );
119
120//------------------------------------------------------------------------------
121
124 index_t
125 comm_splitcount( const Vector< index_t > & aLengths, const proc_t aRoot );
126
127 index_t
128 comm_splitcount( index_t aLength );
129
130//------------------------------------------------------------------------------
131
133 template< typename T >
134 void
135 broadcast( T & aMessage, const proc_t aRoot=0,
136 typename std::enable_if<is_scalar<T>::value>::type* = nullptr )
137 {
138#ifdef BELFEM_MPI
139 if( gComm.size() > 1 )
140 {
141 // Note: std::complex<T> is not supported here; use
142 // send/receive for complex scalar communication.
143 BELFEM_ASSERT( std::is_arithmetic<T>::value,
144 "Can only broadcast arithmetic types." );
145
146 comm_t tCommType = comm_type<T>();
147 comm_check( MPI_Bcast(
148 &aMessage,
149 1,
150 tCommType,
151 aRoot,
152 gComm.world() ) );
153 }
154#endif
155 }
156
158 template< typename T >
159 void
160 broadcast( T * aMessage, const proc_t aRoot, const proc_t aLength )
161 {
162#ifdef BELFEM_MPI
163 if( gComm.size() > 1 )
164 {
165 // Note: std::complex<T> is not supported here; use
166 // send/receive for complex array communication.
167 BELFEM_ASSERT( std::is_arithmetic<T>::value,
168 "Can only broadcast arithmetic types." );
169
170 comm_t tCommType = comm_type<T>();
171
172 comm_check( MPI_Bcast(
173 aMessage,
174 aLength,
175 tCommType,
176 aRoot,
177 gComm.world() ) );
178 }
179#endif
180 }
181
182//------------------------------------------------------------------------------
183
201 template < typename T >
202 void
203 allreduce( const T * aSend, T * aRecv, const int aCount )
204 {
205 // MPI_MAX is meaningless for the complex types comm_type<> also
206 // serves - refuse them at compile time
207 static_assert( std::is_arithmetic< T >::value,
208 "allreduce: MAX-reduction is defined for arithmetic types only" );
209
210 BELFEM_ASSERT( aCount >= 0, "allreduce: negative count %d", aCount );
211
212#ifdef BELFEM_MPI
213 comm_check( MPI_Allreduce(
214 aSend,
215 aRecv,
216 aCount,
217 comm_type< T >(),
218 MPI_MAX,
219 gComm.world() ) );
220#else
221 // serial semantic of a one-rank allreduce: the identity copy
222 if ( aSend != aRecv )
223 {
224 for ( int i = 0; i < aCount; ++i )
225 {
226 aRecv[ i ] = aSend[ i ];
227 }
228 }
229#endif
230 }
231
232//------------------------------------------------------------------------------
233
243 template < typename T >
244 void
245 allreduce_min( const T * aSend, T * aRecv, const int aCount )
246 {
247 static_assert( std::is_arithmetic< T >::value,
248 "allreduce_min: MIN-reduction is defined for arithmetic types only" );
249
250 BELFEM_ASSERT( aCount >= 0, "allreduce_min: negative count %d", aCount );
251
252#ifdef BELFEM_MPI
253 comm_check( MPI_Allreduce(
254 aSend,
255 aRecv,
256 aCount,
257 comm_type< T >(),
258 MPI_MIN,
259 gComm.world() ) );
260#else
261 if ( aSend != aRecv )
262 {
263 for ( int i = 0; i < aCount; ++i )
264 {
265 aRecv[ i ] = aSend[ i ];
266 }
267 }
268#endif
269 }
270
271//==============================================================================
272// SCALARS
273//==============================================================================
274
282 template< typename T >
283 void
284 send( const T aData, const proc_t aTarget=0, typename std::enable_if<is_scalar<T>::value>::type* = nullptr )
285 {
286#ifdef BELFEM_MPI
287 proc_t tMyRank = gComm.rank();
288
289 if( aTarget < comm_size() && tMyRank != aTarget )
290 {
291 MPI_Status tStatus;
292 MPI_Request tRequest;
293
294 comm_check( MPI_Isend( &aData,
295 1,
296 comm_type< T >(),
297 aTarget,
298 comm_tag( aTarget, tMyRank ),
299 gComm.world(),
300 & tRequest ) );
301
302 comm_check( MPI_Wait( &tRequest, &tStatus ) );
303 }
304#endif
305 }
306
307//------------------------------------------------------------------------------
308
313 template< typename T >
314 void
315 receive( T & aData, const proc_t aSource=0, typename std::enable_if<is_scalar<T>::value>::type* = nullptr )
316 {
317#ifdef BELFEM_MPI
318 proc_t tMyRank = gComm.rank();
319
320 proc_t tCommSize = gComm.size();
321
322 if( aSource < tCommSize && tMyRank != aSource )
323 {
324 MPI_Status tStatus;
325 MPI_Request tRequest;
326
327 comm_check( MPI_Irecv( &aData,
328 1,
329 comm_type< T >(),
330 aSource,
331 comm_tag( aSource, tMyRank ),
332 gComm.world(),
333 & tRequest ) );
334
335 comm_check( MPI_Wait( &tRequest, &tStatus ) );
336 }
337#endif
338 }
339
340//==============================================================================
341// RAW ARRAYS
342//==============================================================================
343
346 template< typename T >
347 void
348 send( T * aData, const index_t aLength, const proc_t aTarget )
349 {
350#ifdef BELFEM_MPI
351 proc_t tCommSize = gComm.size();
352
353 proc_t tMyRank = gComm.rank();
354
355 if ( aTarget < tCommSize && aTarget != tMyRank )
356 {
357 MPI_Status tSizeStatus;
358 MPI_Request tSizeRequest;
359
360 int tCommTag = comm_tag( tMyRank, aTarget );
361
362 comm_check( MPI_Isend( &aLength,
363 1,
364 comm_type< index_t >(),
365 aTarget,
366 tCommTag++,
367 gComm.world(),
368 & tSizeRequest ) );
369
370 comm_check( MPI_Wait( &tSizeRequest, &tSizeStatus ) );
371
372 if ( aLength == 0 ) return ;
373
374 Cell< int > tChunkSizes = comm_split( aLength );
375
376 MPI_Status* tStatus = ( MPI_Status* ) malloc( sizeof( MPI_Status ) * tChunkSizes.size() );
377 MPI_Request* tRequest = ( MPI_Request* ) malloc( sizeof( MPI_Request ) * tChunkSizes.size() );
378
379 index_t tCount = 0 ;
380
381 index_t tOffset = 0 ;
382
383 comm_t tCommType = comm_type< T >();
384
385 for ( index_t c : tChunkSizes )
386 {
387 comm_check( MPI_Isend(
388 &aData[ tOffset ],
389 c,
390 tCommType,
391 aTarget,
392 tCommTag,
393 gComm.world(),
394 & tRequest[ tCount++] ) );
395
396 tOffset+= c;
397 }
398
399 comm_check( MPI_Waitall( tCount, tRequest, tStatus ) );
400
401 free( tStatus );
402 free( tRequest );
403 }
404#endif
405 }
406
411 template< typename T >
412 void
413 receive( T * aData, index_t & aLength, const proc_t aSource )
414 {
415#ifdef BELFEM_MPI
416
417 proc_t tCommSize = gComm.size();
418
419 proc_t tMyRank = gComm.rank();
420
421 if ( aSource < tCommSize && aSource != tMyRank )
422 {
423 MPI_Status tSizeStatus;
424 MPI_Request tSizeRequest;
425
426 int tCommTag = comm_tag( aSource, tMyRank );
427
428 index_t tSize = 0;
429
430 comm_check( MPI_Irecv( &tSize,
431 1,
432 comm_type< index_t >(),
433 aSource,
434 tCommTag++,
435 gComm.world(),
436 & tSizeRequest ) );
437
438 comm_check( MPI_Wait( &tSizeRequest, &tSizeStatus ) );
439
440 BELFEM_ERROR( tSize <= aLength, "Datastream too large ( length %u but expect <=%u ).",
441 ( unsigned int ) tSize, ( unsigned int ) aLength );
442
443 aLength = tSize;
444
445 if ( aLength == 0 ) return ;
446
447 Cell< int > tChunkSizes = comm_split( tSize );
448
449 MPI_Status* tStatus = ( MPI_Status* ) malloc( sizeof( MPI_Status ) * tChunkSizes.size() );
450 MPI_Request* tRequest = ( MPI_Request* ) malloc( sizeof( MPI_Request ) * tChunkSizes.size() );
451
452 index_t tCount = 0 ;
453
454 index_t tOffset = 0 ;
455
456 comm_t tCommType = comm_type< T >();
457
458 for ( index_t c : tChunkSizes )
459 {
460 comm_check( MPI_Irecv(
461 &aData[ tOffset ],
462 c,
463 tCommType,
464 aSource,
465 tCommTag,
466 gComm.world(),
467 & tRequest[ tCount++] ) );
468
469 tOffset+= c;
470 }
471
472 comm_check( MPI_Waitall( tCount, tRequest, tStatus ) );
473
474 free( tStatus );
475 free( tRequest );
476 }
477#endif
478 }
479
480//==============================================================================
481// CELLS
482//==============================================================================
483
485 template< typename T >
486 void
487 send( Cell< T > & aData, const proc_t aTarget=0 )
488 {
489#ifdef BELFEM_MPI
490 proc_t tMyRank = gComm.rank();
491
492 if ( tMyRank == aTarget ) return ;
493
494 index_t tSize = aData.size();
495
496 MPI_Status tSizeStatus;
497 MPI_Request tSizeRequest;
498
499 int tCommTag = comm_tag( tMyRank, aTarget );
500
501 comm_check( MPI_Isend( &tSize,
502 1,
503 comm_type< index_t >(),
504 aTarget,
505 tCommTag++,
506 gComm.world(),
507 & tSizeRequest ) );
508
509 comm_check( MPI_Wait( &tSizeRequest, &tSizeStatus ) );
510
511 if ( tSize == 0 ) return ;
512
513 comm_t tCommType = comm_type< T >();
514
515 Cell< int > tChunkSizes = comm_split( aData.size() );
516
517 index_t tOffset = 0 ;
518
519 MPI_Status* tStatus = ( MPI_Status* ) malloc( sizeof( MPI_Status ) * tChunkSizes.size() );
520 MPI_Request* tRequest = ( MPI_Request* ) malloc( sizeof( MPI_Request ) * tChunkSizes.size() );
521
522 index_t tCount = 0 ;
523
524 const T * tData = aData.data();
525
526 for ( index_t c : tChunkSizes )
527 {
528 comm_check( MPI_Isend( &tData[ tOffset ],
529 c,
530 tCommType,
531 aTarget,
532 tCommTag,
533 gComm.world(),
534 &tRequest[ tCount++ ] ) );
535
536 tOffset+= c;
537 }
538
539 comm_check( MPI_Waitall( tCount, tRequest, tStatus ) );
540
541 free( tStatus );
542 free( tRequest );
543#endif
544 }
545
546//------------------------------------------------------------------------------
547
550 template< typename T >
551 void
552 receive( Cell< T > & aData, const proc_t aSource=0 )
553 {
554#ifdef BELFEM_MPI
555
556 proc_t tMyRank = gComm.rank();
557
558 if ( tMyRank == aSource ) return ;
559
560 index_t tSize = 0 ;
561 receive( tSize, aSource );
562
563 aData.set_size( tSize, 0 );
564
565 if ( tSize == 0 ) return ;
566
567 comm_t tCommType = comm_type< T >();
568
569 Cell< int > tChunkSizes = comm_split( aData.size() );
570
571 index_t tOffset = 0 ;
572
573 int tCommTag = comm_tag( aSource, tMyRank ) + 1 ;
574
575 MPI_Status* tStatus = ( MPI_Status* ) malloc( sizeof( MPI_Status ) * tChunkSizes.size() );
576 MPI_Request* tRequest = ( MPI_Request* ) malloc( sizeof( MPI_Request ) * tChunkSizes.size() );
577
578 index_t tCount = 0 ;
579
580 T * tData = aData.data();
581
582 for ( index_t c : tChunkSizes )
583 {
584 comm_check( MPI_Irecv(
585 &tData[ tOffset ],
586 c,
587 tCommType,
588 aSource,
589 tCommTag,
590 gComm.world(),
591 &tRequest[ tCount++ ] ) );
592
593 tOffset+= c;
594 }
595
596 comm_check( MPI_Waitall( tCount, tRequest, tStatus ) );
597
598 free( tStatus );
599 free( tRequest );
600#endif
601 }
602
603//------------------------------------------------------------------------------
604
606 template< typename T >
607 void
608 broadcast( Cell< T > & aData, const proc_t aRoot=0 )
609 {
610#ifdef BELFEM_MPI
611
612 proc_t tMyRank = gComm.rank();
613
614 index_t tSize = tMyRank == aRoot ? aData.size() : 0 ;
615
616 MPI_Request tSizeRequest ;
617 comm_check( MPI_Ibcast(
618 & tSize,
619 1,
620 comm_type< index_t >(),
621 aRoot,
622 gComm.world(),
623 & tSizeRequest ) );
624
625 MPI_Status tSizeStatus ;
626 comm_check( MPI_Wait( &tSizeRequest, &tSizeStatus ) );
627
628 if ( tMyRank != aRoot ) aData.set_size( tSize );
629
630 if ( tSize == 0 ) return ;
631
632 MPI_Request tDataRequest ;
633 comm_check( MPI_Ibcast(
634 aData.data(),
635 tSize,
636 comm_type< T >(),
637 aRoot,
638 gComm.world(),
639 & tDataRequest ) );
640
641 MPI_Status tDataStatus ;
642 comm_check( MPI_Wait( &tDataRequest, &tDataStatus ) );
643
644#endif
645 }
646
647//------------------------------------------------------------------------------
648
655 template< typename T >
656 void
657 distribute( Cell< T > & aData )
658 {
659#ifdef BELFEM_MPI
660 proc_t tMyRank = gComm.rank();
661
662 proc_t tCommSize = gComm.size();
663
664 BELFEM_ASSERT( static_cast< proc_t >( aData.size() )== tCommSize,
665 "Length of cell does not match ( is %u, expect commsize %u ).",
666 ( unsigned int ) aData.size(), ( unsigned int ) tCommSize );
667
668 MPI_Status* tStatus = ( MPI_Status* ) malloc( sizeof( MPI_Status ) * tCommSize );
669 MPI_Request* tRequest = ( MPI_Request* ) malloc( sizeof( MPI_Request ) * tCommSize );
670
671 for( proc_t p=0; p<tCommSize; ++p )
672 {
673 if( p == tMyRank )
674 {
675 tRequest[ p ] = MPI_REQUEST_NULL;
676 continue ;
677 }
678
679 comm_check( MPI_Isend( &aData( p ),
680 1,
681 comm_type< T >(),
682 p,
683 comm_tag( tMyRank, p ),
684 gComm.world(),
685 &tRequest[ p ] ) );
686
687 }
688
689 comm_check( MPI_Waitall( tCommSize, tRequest, tStatus ) );
690
691 free( tStatus );
692 free( tRequest );
693#endif
694 }
695
696//------------------------------------------------------------------------------
697
700 template< typename T >
701 void
702 collect( Cell< T > & aData, const T aMyValue = 0 )
703 {
704#ifdef BELFEM_MPI
705 proc_t tMyRank = gComm.rank();
706
707 proc_t tCommSize = gComm.size();
708
709 aData.set_size( tCommSize, 0 );
710
711 MPI_Status* tStatus = ( MPI_Status* ) malloc( sizeof( MPI_Status ) * tCommSize );
712 MPI_Request* tRequest = ( MPI_Request* ) malloc( sizeof( MPI_Request ) * tCommSize );
713
714 for( proc_t p=0; p<tCommSize; ++p )
715 {
716 if( p == tMyRank )
717 {
718 tRequest[ p ] = MPI_REQUEST_NULL;
719 aData( p ) = aMyValue ;
720 continue ;
721 }
722
723 comm_check( MPI_Irecv( &aData( p ),
724 1,
725 comm_type< T >(),
726 p,
727 comm_tag( tMyRank, p ),
728 gComm.world(),
729 &tRequest[ p ] ) );
730
731 }
732 comm_check( MPI_Waitall( tCommSize, tRequest, tStatus ) );
733
734 free( tStatus );
735 free( tRequest );
736#endif
737 }
738
739//==============================================================================
740// VECTORS
741//==============================================================================
742
744 template< typename T >
745 void
746 send( Vector< T > & aData, const proc_t aTarget=0 )
747 {
748#ifdef BELFEM_MPI
749 proc_t tMyRank = gComm.rank();
750
751 if ( tMyRank == aTarget ) return ;
752
753 index_t tSize = aData.length();
754
755 MPI_Status tSizeStatus;
756 MPI_Request tSizeRequest;
757
758 int tCommTag = comm_tag( tMyRank, aTarget );
759
760 comm_check( MPI_Isend( &tSize,
761 1,
762 comm_type< index_t >(),
763 aTarget,
764 tCommTag++,
765 gComm.world(),
766 & tSizeRequest ) );
767
768 comm_check( MPI_Wait( &tSizeRequest, &tSizeStatus ) );
769
770 if ( tSize == 0 ) return ;
771
772 comm_t tCommType = comm_type< T >();
773
774 Cell< int > tChunkSizes = comm_split( aData.length() );
775
776 index_t tOffset = 0 ;
777
778 MPI_Status* tStatus = ( MPI_Status* ) malloc( sizeof( MPI_Status ) * tChunkSizes.size() );
779 MPI_Request* tRequest = ( MPI_Request* ) malloc( sizeof( MPI_Request ) * tChunkSizes.size() );
780
781 index_t tCount = 0 ;
782
783 const T * tData = aData.data();
784
785 for ( index_t c : tChunkSizes )
786 {
787 comm_check( MPI_Isend( &tData[ tOffset ],
788 c,
789 tCommType,
790 aTarget,
791 tCommTag,
792 gComm.world(),
793 &tRequest[ tCount++ ] ) );
794
795 tOffset+= c;
796 }
797
798 comm_check( MPI_Waitall( tCount, tRequest, tStatus ) );
799
800 free( tStatus );
801 free( tRequest );
802#endif
803 }
804
805//------------------------------------------------------------------------------
806
809 template< typename T >
810 void
811 receive( Vector< T > & aData, const proc_t aSource=0 )
812 {
813#ifdef BELFEM_MPI
814
815 proc_t tMyRank = gComm.rank();
816
817 if ( tMyRank == aSource ) return ;
818
819 index_t tSize = 0 ;
820 receive( tSize, aSource );
821
822 aData.set_size( tSize, 0 );
823
824 if ( tSize == 0 ) return ;
825
826 comm_t tCommType = comm_type< T >();
827
828 Cell< int > tChunkSizes = comm_split( aData.length() );
829
830 index_t tOffset = 0 ;
831
832 int tCommTag = comm_tag( aSource, tMyRank ) + 1 ;
833
834 MPI_Status* tStatus = ( MPI_Status* ) malloc( sizeof( MPI_Status ) * tChunkSizes.size() );
835 MPI_Request* tRequest = ( MPI_Request* ) malloc( sizeof( MPI_Request ) * tChunkSizes.size() );
836
837 index_t tCount = 0 ;
838
839 T * tData = aData.data();
840
841 for ( index_t c : tChunkSizes )
842 {
843 comm_check( MPI_Irecv(
844 &tData[ tOffset ],
845 c,
846 tCommType,
847 aSource,
848 tCommTag,
849 gComm.world(),
850 &tRequest[ tCount++ ] ) );
851
852 tOffset+= c;
853 }
854
855 comm_check( MPI_Waitall( tCount, tRequest, tStatus ) );
856
857 free( tStatus );
858 free( tRequest );
859#endif
860 }
861
862//------------------------------------------------------------------------------
863
865 template< typename T >
866 void
867 broadcast( Vector< T > & aData, const proc_t aRoot=0 )
868 {
869#ifdef BELFEM_MPI
870
871 proc_t tMyRank = gComm.rank();
872
873 index_t tSize = tMyRank == aRoot ? aData.length() : 0 ;
874
875 MPI_Request tSizeRequest ;
876 comm_check( MPI_Ibcast(
877 & tSize,
878 1,
879 comm_type< index_t >(),
880 aRoot,
881 gComm.world(),
882 & tSizeRequest ) );
883
884 MPI_Status tSizeStatus ;
885 comm_check( MPI_Wait( &tSizeRequest, &tSizeStatus ) );
886
887 if ( tMyRank != aRoot ) aData.set_size( tSize );
888
889 if ( tSize == 0 ) return ;
890
891 MPI_Request tDataRequest ;
892 comm_check( MPI_Ibcast(
893 aData.data(),
894 tSize,
895 comm_type< T >(),
896 aRoot,
897 gComm.world(),
898 & tDataRequest ) );
899
900 MPI_Status tDataStatus ;
901 comm_check( MPI_Wait( &tDataRequest, &tDataStatus ) );
902
903#endif
904 }
905
906//------------------------------------------------------------------------------
907
909 template< typename T >
910 void
911 distribute( Vector< T > & aData )
912 {
913#ifdef BELFEM_MPI
914 proc_t tMyRank = gComm.rank();
915
916 proc_t tCommSize = gComm.size();
917
918 BELFEM_ASSERT( static_cast< proc_t>( aData.length() )== tCommSize,
919 "Length of vector does not match ( is %u, expect commsize %u ).",
920 ( unsigned int ) aData.length(), ( unsigned int ) tCommSize );
921
922 MPI_Status* tStatus = ( MPI_Status* ) malloc( sizeof( MPI_Status ) * tCommSize );
923 MPI_Request* tRequest = ( MPI_Request* ) malloc( sizeof( MPI_Request ) * tCommSize );
924
925 for( proc_t p=0; p<tCommSize; ++p )
926 {
927 if( p == tMyRank )
928 {
929 tRequest[ p ] = MPI_REQUEST_NULL;
930 continue ;
931 }
932
933 comm_check( MPI_Isend( &aData( p ),
934 1,
935 comm_type< T >(),
936 p,
937 comm_tag( tMyRank, p ),
938 gComm.world(),
939 &tRequest[ p ] ) );
940
941 }
942
943 comm_check( MPI_Waitall( tCommSize, tRequest, tStatus ) );
944
945 free( tStatus );
946 free( tRequest );
947#endif
948 }
949
950//------------------------------------------------------------------------------
951
954 template< typename T >
955 void
956 collect( Vector< T > & aData, const T aMyValue = 0 )
957 {
958#ifdef BELFEM_MPI
959 proc_t tMyRank = gComm.rank();
960
961 proc_t tCommSize = gComm.size();
962
963 aData.set_size( tCommSize, 0 );
964
965 MPI_Status* tStatus = ( MPI_Status* ) malloc( sizeof( MPI_Status ) * tCommSize );
966 MPI_Request* tRequest = ( MPI_Request* ) malloc( sizeof( MPI_Request ) * tCommSize );
967
968 for( proc_t p=0; p<tCommSize; ++p )
969 {
970 if( p == tMyRank )
971 {
972 tRequest[ p ] = MPI_REQUEST_NULL;
973 aData( p ) = aMyValue ;
974 continue ;
975 }
976
977 comm_check( MPI_Irecv( &aData( p ),
978 1,
979 comm_type< T >(),
980 p,
981 comm_tag( tMyRank, p ),
982 gComm.world(),
983 &tRequest[ p ] ) );
984
985 }
986 comm_check( MPI_Waitall( tCommSize, tRequest, tStatus ) );
987
988 free( tStatus );
989 free( tRequest );
990#endif
991 }
992
993//------------------------------------------------------------------------------
994
996 template< typename T >
997 void
998 distribute( Cell< Vector< T > > & aData )
999 {
1000#ifdef BELFEM_MPI
1001
1002 proc_t tMyRank = gComm.rank();
1003
1004 proc_t tCommSize = gComm.size();
1005
1006 BELFEM_ASSERT( static_cast< proc_t>( aData.size() )== tCommSize,
1007 "Length of data container does not match ( is %u, expect commsize %u ).",
1008 ( unsigned int ) aData.size(), ( unsigned int ) tCommSize );
1009
1010 Vector< index_t > tSizes( tCommSize, 0 );
1011 for ( proc_t p=0; p<tCommSize; ++p )
1012 {
1013 tSizes( p ) = aData( p ).length();
1014 }
1015 index_t tCount = comm_splitcount( tSizes, tMyRank );
1016
1017 distribute( tSizes );
1018
1019 MPI_Request* tRequest = ( MPI_Request * ) malloc( sizeof( MPI_Request ) * tCount );
1020 MPI_Status* tStatus = ( MPI_Status * ) malloc( sizeof( MPI_Status ) * tCount );
1021
1022 tCount = 0 ;
1023
1024 comm_t tCommType = comm_type< T >();
1025
1026 for( proc_t p=0; p<tCommSize; ++p )
1027 {
1028 if ( p == tMyRank || aData( p ).length() == 0 ) continue ;
1029
1030 const T * tData = aData( p ).data();
1031
1032 Cell< int > tChunkSizes = comm_split( aData( p ).length() );
1033
1034 index_t tOffset = 0 ;
1035
1036 int tCommTag = comm_tag( tMyRank, p ) + 1 ;
1037
1038 for ( index_t c : tChunkSizes )
1039 {
1040 comm_check( MPI_Isend( &tData[ tOffset ],
1041 c,
1042 tCommType,
1043 p,
1044 tCommTag,
1045 gComm.world(),
1046 &tRequest[ tCount++ ] ) );
1047
1048 tOffset+= c;
1049 }
1050 }
1051
1052 comm_check( MPI_Waitall( tCount, tRequest, tStatus ) );
1053
1054 free( tStatus );
1055 free( tRequest );
1056
1057#endif
1058 }
1059
1061 template< typename T >
1062 void
1063 distribute( Cell< Cell< T > > & aData )
1064 {
1065#ifdef BELFEM_MPI
1066
1067 proc_t tMyRank = gComm.rank();
1068
1069 proc_t tCommSize = gComm.size();
1070
1071 BELFEM_ASSERT( static_cast< proc_t>( aData.size() )== tCommSize,
1072 "Length of data container does not match ( is %u, expect commsize %u ).",
1073 ( unsigned int ) aData.size(), ( unsigned int ) tCommSize );
1074
1075 Vector< index_t > tSizes( tCommSize, 0 );
1076 for ( proc_t p=0; p<tCommSize; ++p )
1077 {
1078 tSizes( p ) = aData( p ).size();
1079 }
1080 index_t tCount = comm_splitcount( tSizes, tMyRank );
1081
1082 distribute( tSizes );
1083
1084 MPI_Request* tRequest = ( MPI_Request * ) malloc( sizeof( MPI_Request ) * tCount );
1085 MPI_Status* tStatus = ( MPI_Status * ) malloc( sizeof( MPI_Status ) * tCount );
1086
1087 tCount = 0 ;
1088
1089 comm_t tCommType = comm_type< T >();
1090
1091 for( proc_t p=0; p<tCommSize; ++p )
1092 {
1093 if ( p == tMyRank || aData( p ).size() == 0 ) continue ;
1094
1095 const T * tData = aData( p ).data();
1096
1097 Cell< int > tChunkSizes = comm_split( aData( p ).size() );
1098
1099 index_t tOffset = 0 ;
1100
1101 int tCommTag = comm_tag( tMyRank, p ) + 1 ;
1102
1103 for ( index_t c : tChunkSizes )
1104 {
1105 comm_check( MPI_Isend( &tData[ tOffset ],
1106 c,
1107 tCommType,
1108 p,
1109 tCommTag,
1110 gComm.world(),
1111 &tRequest[ tCount++ ] ) );
1112
1113 tOffset+= c;
1114 }
1115 }
1116
1117 comm_check( MPI_Waitall( tCount, tRequest, tStatus ) );
1118
1119 free( tStatus );
1120 free( tRequest );
1121
1122#endif
1123 }
1124
1125//------------------------------------------------------------------------------
1126
1130 template< typename T, typename U >
1131 void
1132 distribute( const T * aData, const Vector< U > & aOffsets )
1133 {
1134#ifdef BELFEM_MPI
1135
1136 proc_t tMyRank = gComm.rank();
1137
1138 proc_t tCommSize = gComm.size();
1139
1140 BELFEM_ASSERT( static_cast< proc_t>( aOffsets.length() ) == tCommSize + 1,
1141 "Length of offset container does not match ( is %u, expect commsize %u ).",
1142 ( unsigned int ) aOffsets.length(), ( unsigned int ) tCommSize + 1 );
1143
1144 Vector< index_t > tSizes( tCommSize, 0 );
1145
1146 for ( proc_t p=0; p<tCommSize; ++p )
1147 {
1148 tSizes( p ) = aOffsets( p+1 ) - aOffsets( p );
1149 }
1150 index_t tCount = comm_splitcount( tSizes, tMyRank );
1151
1152 distribute( tSizes );
1153
1154 MPI_Request* tRequest = ( MPI_Request * ) malloc( sizeof( MPI_Request ) * tCount );
1155 MPI_Status* tStatus = ( MPI_Status * ) malloc( sizeof( MPI_Status ) * tCount );
1156
1157 tCount = 0 ;
1158
1159 comm_t tCommType = comm_type< typename std::remove_const<T>::type >();
1160
1161 for( proc_t p=0; p<tCommSize; ++p )
1162 {
1163 if ( p == tMyRank || tSizes( p ) == 0 ) continue ;
1164
1165 Cell< int > tChunkSizes = comm_split( tSizes( p ) );
1166
1167 index_t tOffset = aOffsets( p ) ;
1168
1169 int tCommTag = comm_tag( tMyRank, p ) + 1 ;
1170
1171 for ( index_t c : tChunkSizes )
1172 {
1173 comm_check( MPI_Isend( aData + tOffset,
1174 c,
1175 tCommType,
1176 p,
1177 tCommTag,
1178 gComm.world(),
1179 &tRequest[ tCount++ ] ) );
1180
1181 tOffset+= c;
1182 }
1183 }
1184
1185 comm_check( MPI_Waitall( tCount, tRequest, tStatus ) );
1186
1187 free( tStatus );
1188 free( tRequest );
1189
1190#endif
1191 }
1192
1195 template< typename T, typename U >
1196 void
1197 collect( T * aData, const Vector< U > & aOffsets )
1198 {
1199#ifdef BELFEM_MPI
1200 proc_t tMyRank = gComm.rank();
1201
1202 proc_t tCommSize = gComm.size();
1203
1204 Vector< index_t > tSizes ;
1205 collect( tSizes );
1206 // on rank 0, slot 0 is this rank's own block and is already in place;
1207 // the loops below start at p = 1 and never read tSizes( 0 )
1208 tSizes( 0 ) = aOffsets( 0 );
1209
1210 index_t tCount = 0 ;
1211
1212 for ( proc_t p=1; p<tCommSize; ++p )
1213 {
1214 tCount += comm_split( tSizes( p ) ).size();
1215 }
1216
1217 MPI_Status* tStatus = ( MPI_Status* ) malloc( sizeof( MPI_Status ) * tCount );
1218 MPI_Request* tRequest = ( MPI_Request* ) malloc( sizeof( MPI_Request ) * tCount );
1219
1220 tCount = 0 ;
1221
1222 comm_t tCommType = comm_type< T >();
1223
1224 for( proc_t p=1; p<tCommSize; ++p )
1225 {
1226 index_t tSize = tSizes( p );
1227
1228 if ( tSize == 0 ) continue;
1229
1230 Cell< int > tChunkSizes = comm_split( tSize );
1231
1232 index_t tOffset = aOffsets( p ) ;
1233
1234 int tCommTag = comm_tag( tMyRank, p ) + 1 ;
1235
1236 for ( index_t c : tChunkSizes )
1237 {
1238 comm_check( MPI_Irecv( &aData[ tOffset ],
1239 c,
1240 tCommType,
1241 p,
1242 tCommTag,
1243 gComm.world(),
1244 &tRequest[ tCount++ ] ) );
1245
1246 tOffset+= c;
1247 }
1248 }
1249
1250 comm_check( MPI_Waitall( tCount, tRequest, tStatus ) );
1251
1252 free( tStatus );
1253 free( tRequest );
1254#endif
1255 }
1256//------------------------------------------------------------------------------
1257
1260 template< typename T >
1261 void
1262 collect( Cell< Vector< T > > & aData, Vector< T > aMyData={} )
1263 {
1264#ifdef BELFEM_MPI
1265
1266 proc_t tMyRank = gComm.rank();
1267
1268 proc_t tCommSize = gComm.size();
1269
1270 aData.set_size( tCommSize, {} );
1271
1272 Vector< index_t > tSizes ;
1273 collect( tSizes );
1274
1275 index_t tCount = 0 ;
1276 for ( proc_t p=0; p<tCommSize; ++p )
1277 {
1278 if ( p == tMyRank )
1279 {
1280 if( aMyData.length() > 0 )
1281 {
1282 aData( p ).vector_data() = std::move( aMyData.vector_data() ) ;
1283 }
1284 continue ;
1285 }
1286
1287 aData( p ).set_size( tSizes( p ), 0 );
1288
1289 tCount += comm_split( tSizes( p ) ).size();
1290 }
1291
1292 MPI_Status* tStatus = ( MPI_Status* ) malloc( sizeof( MPI_Status ) * tCount );
1293 MPI_Request* tRequest = ( MPI_Request* ) malloc( sizeof( MPI_Request ) * tCount );
1294
1295 tCount = 0 ;
1296
1297 comm_t tCommType = comm_type< T >();
1298
1299 for( proc_t p=0; p<tCommSize; ++p )
1300 {
1301 index_t tSize = tSizes( p );
1302
1303 if ( p == tMyRank || tSize == 0 ) continue;
1304
1305 T * tData = aData( p ).data() ;
1306
1307 Cell< int > tChunkSizes = comm_split( tSize );
1308
1309 index_t tOffset = 0 ;
1310
1311 int tCommTag = comm_tag( tMyRank, p ) + 1 ;
1312
1313 for ( index_t c : tChunkSizes )
1314 {
1315 comm_check( MPI_Irecv( &tData[ tOffset ],
1316 c,
1317 tCommType,
1318 p,
1319 tCommTag,
1320 gComm.world(),
1321 &tRequest[ tCount++ ] ) );
1322
1323 tOffset+= c;
1324 }
1325 }
1326
1327 comm_check( MPI_Waitall( tCount, tRequest, tStatus ) );
1328
1329 free( tStatus );
1330 free( tRequest );
1331
1332#endif
1333 }
1334
1335//------------------------------------------------------------------------------
1336
1342 template< typename T >
1343 void
1344 share( Vector< T > & aData )
1345 {
1346#ifdef BELFEM_MPI
1347 proc_t tMyRank = gComm.rank();
1348
1349 proc_t tCommSize = gComm.size();
1350
1351 index_t tSize = aData.length();
1352
1353 MPI_Request* tRequest = ( MPI_Request * ) malloc( sizeof( MPI_Request ) * ( tCommSize-1 ) );
1354 MPI_Status* tStatus = ( MPI_Status * ) malloc( sizeof( MPI_Status ) * ( tCommSize-1 ) );
1355
1356 comm_t tIndex_t = comm_type< index_t >();
1357
1358 comm_t tCommType = comm_type< T >();
1359
1360 index_t tCount = 0 ;
1361
1362 for ( proc_t p=0; p<tCommSize; ++p )
1363 {
1364 if ( p == tMyRank ) continue;
1365
1366 int tCommTag = comm_tag( tMyRank, p );
1367
1368 comm_check( MPI_Isend( &tSize,
1369 1,
1370 tIndex_t,
1371 p,
1372 tCommTag,
1373 gComm.world(),
1374 &tRequest[ tCount ++ ] ) );
1375 }
1376
1377 comm_check( MPI_Waitall( tCount, tRequest, tStatus ) );
1378
1379 free( tStatus );
1380 free( tRequest );
1381
1382 tCount = comm_splitcount( tSize );
1383
1384 tRequest = ( MPI_Request * ) malloc( sizeof( MPI_Request ) * tCount );
1385 tStatus = ( MPI_Status * ) malloc( sizeof( MPI_Status ) * tCount );
1386
1387 if ( aData.length() > 0 )
1388 {
1389 tCount = 0 ;
1390 for( proc_t p=0; p<tCommSize; ++p )
1391 {
1392 if ( p == tMyRank ) continue ;
1393
1394 const T * tData = aData.data();
1395
1396 Cell< int > tChunkSizes = comm_split( aData.length() );
1397
1398 index_t tOffset = 0 ;
1399
1400 int tCommTag = comm_tag( tMyRank, p ) + 1 ;
1401
1402 for ( index_t c : tChunkSizes )
1403 {
1404 comm_check( MPI_Isend( &tData[ tOffset ],
1405 c,
1406 tCommType,
1407 p,
1408 tCommTag,
1409 gComm.world(),
1410 &tRequest[ tCount++ ] ) );
1411
1412 tOffset+= c;
1413 }
1414 }
1415
1416 comm_check( MPI_Waitall( tCount, tRequest, tStatus ) );
1417 }
1418
1419 free( tStatus );
1420 free( tRequest );
1421#endif
1422 }
1423
1424//------------------------------------------------------------------------------
1425
1427 template< typename T >
1428 void
1429 share( Cell< T > & aData )
1430 {
1431#ifdef BELFEM_MPI
1432 proc_t tMyRank = gComm.rank();
1433
1434 proc_t tCommSize = gComm.size();
1435
1436 index_t tSize = aData.size();
1437
1438 MPI_Request* tRequest = ( MPI_Request * ) malloc( sizeof( MPI_Request ) * ( tCommSize-1 ) );
1439 MPI_Status* tStatus = ( MPI_Status * ) malloc( sizeof( MPI_Status ) * ( tCommSize-1 ) );
1440
1441 comm_t tIndex_t = comm_type< index_t >();
1442
1443 comm_t tCommType = comm_type< T >();
1444
1445 index_t tCount = 0 ;
1446
1447 for ( proc_t p=0; p<tCommSize; ++p )
1448 {
1449 if ( p == tMyRank ) continue;
1450
1451 int tCommTag = comm_tag( tMyRank, p );
1452
1453 comm_check( MPI_Isend( &tSize,
1454 1,
1455 tIndex_t,
1456 p,
1457 tCommTag,
1458 gComm.world(),
1459 &tRequest[ tCount ++ ] ) );
1460 }
1461
1462 comm_check( MPI_Waitall( tCount, tRequest, tStatus ) );
1463
1464 free( tStatus );
1465 free( tRequest );
1466
1467 tCount = comm_splitcount( tSize );
1468
1469 tRequest = ( MPI_Request * ) malloc( sizeof( MPI_Request ) * tCount );
1470 tStatus = ( MPI_Status * ) malloc( sizeof( MPI_Status ) * tCount );
1471
1472 if ( aData.size() > 0 )
1473 {
1474 tCount = 0 ;
1475 for( proc_t p=0; p<tCommSize; ++p )
1476 {
1477 if ( p == tMyRank ) continue ;
1478
1479 const T * tData = aData.data();
1480
1481 Cell< int > tChunkSizes = comm_split( aData.size() );
1482
1483 index_t tOffset = 0 ;
1484
1485 int tCommTag = comm_tag( tMyRank, p ) + 1 ;
1486
1487 for ( index_t c : tChunkSizes )
1488 {
1489 comm_check( MPI_Isend( &tData[ tOffset ],
1490 c,
1491 tCommType,
1492 p,
1493 tCommTag,
1494 gComm.world(),
1495 &tRequest[ tCount++ ] ) );
1496
1497 tOffset+= c;
1498 }
1499 }
1500
1501 comm_check( MPI_Waitall( tCount, tRequest, tStatus ) );
1502 }
1503
1504 free( tStatus );
1505 free( tRequest );
1506#endif
1507 }
1508
1509//==============================================================================
1510// MATRICES
1511//==============================================================================
1512
1514 template< typename T >
1515 void
1516 broadcast( Matrix< T > & aData, const proc_t aRoot=0 )
1517 {
1518#ifdef BELFEM_MPI
1519
1520 proc_t tMyRank = gComm.rank();
1521
1522 // We transmit spacing()*n_cols rather than n_rows*n_cols so that the
1523 // raw contiguous buffer (including backend padding, e.g. Blaze
1524 // SIMD alignment) is transferred as-is. This avoids disassembling
1525 // and reassembling the matrix. Both ranks use the same binary, so
1526 // padding layout is identical for a given (rows, cols) pair.
1527 // Never use capacity() here: a matrix that shrank keeps its old,
1528 // larger allocation, and transmitting that count overflows the
1529 // exact-fit buffer on the receiving side.
1530 index_t tSize[ 3 ];
1531
1532 tSize[ 0 ] = aData.n_rows();
1533 tSize[ 1 ] = aData.n_cols();
1534 tSize[ 2 ] = aData.spacing() * aData.n_cols();
1535
1536 MPI_Request tSizeRequest ;
1537 comm_check( MPI_Ibcast(
1538 & tSize,
1539 3,
1540 comm_type< index_t >(),
1541 aRoot,
1542 gComm.world(),
1543 & tSizeRequest ) );
1544
1545 MPI_Status tSizeStatus ;
1546 comm_check( MPI_Wait( &tSizeRequest, &tSizeStatus ) );
1547
1548 if ( tMyRank != aRoot ) aData.set_size( tSize[ 0 ], tSize[ 1 ] );
1549
1550 if( tSize[ 0 ] == 0 || tSize[ 1 ] == 0 ) return ;
1551
1552 BELFEM_ERROR( aData.capacity() >= tSize[ 2 ],
1553 "broadcast( Matrix ) : transfer length %lu exceeds local buffer capacity %lu",
1554 ( long unsigned int ) tSize[ 2 ],
1555 ( long unsigned int ) aData.capacity() );
1556
1557 // broadcast is unchunked: the count must fit into MPI's int
1558 BELFEM_ERROR( tSize[ 2 ] <= ( index_t ) std::numeric_limits< int >::max(),
1559 "broadcast( Matrix ) : matrix too large for unchunked broadcast, use send/receive" );
1560
1561 MPI_Request tDataRequest ;
1562 comm_check( MPI_Ibcast(
1563 aData.data(),
1564 tSize[ 2 ],
1565 comm_type< T >(),
1566 aRoot,
1567 gComm.world(),
1568 & tDataRequest ) );
1569
1570 MPI_Status tDataStatus ;
1571 comm_check( MPI_Wait( &tDataRequest, &tDataStatus ) );
1572
1573#endif
1574 }
1575
1577 template< typename T >
1578 void
1579 send( Matrix< T > & aData, const proc_t aTarget=0 )
1580 {
1581#ifdef BELFEM_MPI
1582
1583 proc_t tMyRank = gComm.rank();
1584
1585 if ( tMyRank == aTarget ) return ;
1586
1587 // The transfer length is spacing() * n_cols(), not capacity().
1588 // See broadcast( Matrix ) for the padding-layout rationale.
1589 index_t tSize[ 3 ];
1590
1591 tSize[ 0 ] = aData.n_rows();
1592 tSize[ 1 ] = aData.n_cols();
1593 tSize[ 2 ] = aData.spacing() * aData.n_cols();
1594
1595 MPI_Status tSizeStatus;
1596 MPI_Request tSizeRequest;
1597
1598 int tCommTag = comm_tag( tMyRank, aTarget );
1599
1600 comm_check( MPI_Isend( &tSize,
1601 3,
1602 comm_type< index_t >(),
1603 aTarget,
1604 tCommTag++,
1605 gComm.world(),
1606 & tSizeRequest ) );
1607
1608 comm_check( MPI_Wait( &tSizeRequest, &tSizeStatus ) );
1609
1610 if( tSize[ 0 ] == 0 || tSize[ 1 ] == 0 ) return ;
1611
1612 comm_t tCommType = comm_type< T >();
1613
1614 Cell< int > tChunkSizes = comm_split( tSize[ 2 ] );
1615
1616 index_t tOffset = 0 ;
1617
1618 MPI_Status* tStatus = ( MPI_Status* ) malloc( sizeof( MPI_Status ) * tChunkSizes.size() );
1619 MPI_Request* tRequest = ( MPI_Request* ) malloc( sizeof( MPI_Request ) * tChunkSizes.size() );
1620
1621 index_t tCount = 0 ;
1622
1623 const T * tData = aData.data();
1624
1625 for ( index_t c : tChunkSizes )
1626 {
1627 comm_check( MPI_Isend( &tData[ tOffset ],
1628 c,
1629 tCommType,
1630 aTarget,
1631 tCommTag,
1632 gComm.world(),
1633 &tRequest[ tCount++ ] ) );
1634
1635 tOffset+= c;
1636 }
1637
1638 comm_check( MPI_Waitall( tCount, tRequest, tStatus ) );
1639
1640 free( tStatus );
1641 free( tRequest );
1642#endif
1643 }
1644
1645//------------------------------------------------------------------------------
1646
1650 template< typename T >
1651 void
1652 receive( Matrix< T > & aData, const proc_t aSource=0 )
1653 {
1654#ifdef BELFEM_MPI
1655
1656 proc_t tMyRank = gComm.rank();
1657
1658 if ( tMyRank == aSource ) return ;
1659
1660 MPI_Status tSizeStatus;
1661 MPI_Request tSizeRequest;
1662
1663 int tCommTag = comm_tag( tMyRank, aSource );
1664
1665 // tSize[2] is the sender's spacing()*n_cols, not n_rows*n_cols. We
1666 // transfer the raw contiguous buffer including backend padding
1667 // (e.g. Blaze SIMD alignment) as-is, avoiding matrix disassembly/
1668 // reassembly. Both ranks use the same binary, so padding layout is
1669 // identical for a given (rows, cols) pair, and after set_size below
1670 // our own buffer holds at least spacing()*n_cols elements.
1671 index_t tSize[ 3 ];
1672
1673 comm_check( MPI_Irecv( tSize,
1674 3,
1675 comm_type< index_t >(),
1676 aSource,
1677 tCommTag++,
1678 gComm.world(),
1679 & tSizeRequest ) );
1680
1681 comm_check( MPI_Wait( &tSizeRequest, &tSizeStatus ) );
1682
1683 aData.set_size( tSize[ 0 ], tSize[ 1 ] );
1684
1685 if ( tSize[ 0 ] == 0 || tSize[ 1 ] == 0 ) return ;
1686
1687 BELFEM_ERROR( aData.capacity() >= tSize[ 2 ],
1688 "receive( Matrix ) : transfer length %lu exceeds local buffer capacity %lu",
1689 ( long unsigned int ) tSize[ 2 ],
1690 ( long unsigned int ) aData.capacity() );
1691
1692 comm_t tCommType = comm_type< T >();
1693
1694 Cell< int > tChunkSizes = comm_split( tSize[ 2 ] );
1695
1696 index_t tOffset = 0 ;
1697
1698 MPI_Request* tRequest = ( MPI_Request* ) malloc( sizeof( MPI_Request ) * tChunkSizes.size() );
1699 MPI_Status* tStatus = ( MPI_Status* ) malloc( sizeof( MPI_Status ) * tChunkSizes.size() );
1700
1701 index_t tCount = 0 ;
1702
1703 T * tData = aData.data();
1704
1705 for ( index_t c : tChunkSizes )
1706 {
1707 comm_check( MPI_Irecv( &tData[ tOffset ],
1708 c,
1709 tCommType,
1710 aSource,
1711 tCommTag,
1712 gComm.world(),
1713 &tRequest[ tCount++ ] ) );
1714
1715 tOffset+= c;
1716 }
1717
1718 comm_check( MPI_Waitall( tCount, tRequest, tStatus ) );
1719
1720 free( tStatus );
1721 free( tRequest );
1722#endif
1723 }
1724
1725//------------------------------------------------------------------------------
1726
1728 template< typename T >
1729 void
1730 distribute( Cell< Matrix< T > > & aData )
1731 {
1732#ifdef BELFEM_MPI
1733
1734 proc_t tMyRank = gComm.rank();
1735
1736 proc_t tCommSize = gComm.size();
1737
1738 BELFEM_ASSERT( static_cast< proc_t>( aData.size() )== tCommSize,
1739 "Length of data container does not match ( is %u, expect commsize %u ).",
1740 ( unsigned int ) aData.size(), ( unsigned int ) tCommSize );
1741
1742 index_t * tSizes = ( index_t * ) malloc( sizeof( index_t ) * tCommSize * 3 );
1743 index_t tCount = 0 ;
1744 for ( proc_t p=0; p<tCommSize; ++p )
1745 {
1746 tSizes[ tCount++ ] = aData( p ).n_rows();
1747 tSizes[ tCount++ ] = aData( p ).n_cols();
1748
1749 // transfer length: padded footprint of the current shape,
1750 // NOT capacity() ( which can be stale-large after a shrink )
1751 tSizes[ tCount++ ] = aData( p ).spacing() * aData( p ).n_cols();
1752 }
1753
1754 MPI_Request* tSizeRequest = ( MPI_Request* ) malloc( sizeof( MPI_Request ) * tCommSize );
1755 MPI_Status* tSizeStatus = ( MPI_Status* ) malloc( sizeof( MPI_Status ) * tCommSize );
1756
1757 index_t tOffset = 0 ;
1758 for ( proc_t p=0; p<tCommSize; ++p )
1759 {
1760 if ( p == tMyRank )
1761 {
1762 tSizeRequest[ p ] = MPI_REQUEST_NULL;
1763 tOffset += 3 ;
1764 continue;
1765 }
1766
1767 comm_check( MPI_Isend( &tSizes[ tOffset ],
1768 3,
1769 comm_type< index_t >(),
1770 p,
1771 comm_tag( tMyRank, p ),
1772 gComm.world(),
1773 &tSizeRequest[ p ] ) );
1774
1775 tOffset += 3 ;
1776 }
1777
1778 comm_check( MPI_Waitall( tCommSize, tSizeRequest, tSizeStatus ) );
1779 free( tSizeRequest );
1780 free( tSizeStatus );
1781
1782 tOffset = 2 ;
1783 tCount = 0 ;
1784 for ( proc_t p=0; p<tCommSize; ++p )
1785 {
1786 if ( p == tMyRank || aData( p ).n_rows() == 0 || aData( p ).n_cols() == 0 )
1787 {
1788 tOffset += 3 ;
1789 continue ;
1790 }
1791 tCount += comm_split( tSizes[ tOffset ] ).size();
1792 tOffset += 3 ;
1793 }
1794
1795 free( tSizes );
1796
1797 MPI_Status* tStatus = ( MPI_Status * ) malloc( sizeof( MPI_Status ) * tCount );
1798 MPI_Request* tRequest = ( MPI_Request * ) malloc( sizeof( MPI_Request ) * tCount );
1799
1800 tCount = 0 ;
1801
1802 comm_t tCommType = comm_type< T >();
1803
1804 for( proc_t p=0; p<tCommSize; ++p )
1805 {
1806 if ( p == tMyRank || aData( p ).n_rows() == 0 || aData( p ).n_cols() == 0 )
1807 {
1808 continue;
1809 }
1810
1811 const T * tData = aData( p ).data();
1812
1813 // compute the chunks for this message ( same transfer length
1814 // as announced in tSizes above — never capacity() )
1815 Cell< int > tChunkSizes = comm_split( aData( p ).spacing() * aData( p ).n_cols() );
1816
1817 index_t tDataOffset = 0 ;
1818 int tCommTag = comm_tag( tMyRank, p ) + 1 ;
1819
1820 for ( index_t c : tChunkSizes )
1821 {
1822 comm_check( MPI_Isend(
1823 &tData[ tDataOffset ],
1824 c,
1825 tCommType,
1826 p,
1827 tCommTag,
1828 gComm.world(),
1829 &tRequest[ tCount++ ] ) );
1830
1831 tDataOffset+= c;
1832 }
1833 }
1834
1835 comm_check( MPI_Waitall( tCount, tRequest, tStatus ) );
1836
1837 free( tStatus );
1838 free( tRequest );
1839
1840#endif
1841 }
1842
1843//------------------------------------------------------------------------------
1844
1847 template< typename T >
1848 void
1849 collect( Cell< Matrix< T > > & aData )
1850 {
1851#ifdef BELFEM_MPI
1852
1853 proc_t tMyRank = gComm.rank();
1854
1855 proc_t tCommSize = gComm.size();
1856
1857 aData.set_size( tCommSize, {} );
1858
1859 index_t * tSizes = ( index_t * ) malloc( sizeof( index_t ) * tCommSize * 3 );
1860
1861 MPI_Request* tSizeRequest = ( MPI_Request* ) malloc( sizeof( MPI_Request ) * tCommSize );
1862 MPI_Status* tSizeStatus = ( MPI_Status* ) malloc( sizeof( MPI_Status ) * tCommSize );
1863
1864 index_t tOffset = 0 ;
1865
1866 for ( proc_t p=0; p<tCommSize; ++p )
1867 {
1868 if ( p == tMyRank )
1869 {
1870 tSizeRequest[ p ] = MPI_REQUEST_NULL;
1871 tSizes[ tOffset++ ] = 0 ;
1872 tSizes[ tOffset++ ] = 0 ;
1873 tSizes[ tOffset++ ] = 0 ;
1874 continue ;
1875 }
1876
1877 comm_check( MPI_Irecv( &tSizes[ tOffset ],
1878 3,
1879 comm_type< index_t >(),
1880 p,
1881 comm_tag( tMyRank, p ),
1882 gComm.world(),
1883 &tSizeRequest[ p ] ) );
1884
1885 tOffset += 3 ;
1886 }
1887
1888 comm_check( MPI_Waitall( tCommSize, tSizeRequest, tSizeStatus ) );
1889
1890 free( tSizeStatus );
1891 free( tSizeRequest );
1892
1893 index_t tCount = 0 ;
1894 tOffset = 0 ;
1895
1896 for ( proc_t p=0; p<tCommSize; ++p )
1897 {
1898 if ( p == tMyRank )
1899 {
1900 tOffset += 3 ;
1901 continue ;
1902 }
1903
1904 aData( p ).set_size( tSizes[ tOffset ], tSizes[ tOffset+1 ] );
1905
1906 if ( tSizes[ tOffset ] == 0 || tSizes[ tOffset+1 ] == 0 )
1907 {
1908 tOffset += 3 ;
1909 continue ;
1910 }
1911
1912 tCount += comm_split( tSizes[ tOffset + 2 ] ).size();
1913 tOffset += 3 ;
1914 }
1915
1916 MPI_Status* tStatus = ( MPI_Status* ) malloc( sizeof( MPI_Status ) * tCount );
1917 MPI_Request* tRequest = ( MPI_Request* ) malloc( sizeof( MPI_Request ) * tCount );
1918
1919 comm_t tCommType = comm_type< T >();
1920
1921 tCount = 0 ;
1922
1923 tOffset = 0 ;
1924
1925 for( proc_t p=0; p<tCommSize; ++p )
1926 {
1927 if ( p == tMyRank )
1928 {
1929 tOffset += 3 ;
1930 continue;
1931 }
1932 if ( aData( p ).n_rows() == 0 || aData( p ).n_cols() == 0 )
1933 {
1934 tOffset += 3 ;
1935 continue ;
1936 }
1937
1938 BELFEM_ERROR( aData( p ).capacity() >= tSizes[ tOffset + 2 ],
1939 "collect( Matrix ) : transfer length %lu from proc %u exceeds local buffer capacity %lu",
1940 ( long unsigned int ) tSizes[ tOffset + 2 ],
1941 ( unsigned int ) p,
1942 ( long unsigned int ) aData( p ).capacity() );
1943
1944 Cell< int > tChunkSizes = comm_split( tSizes[ tOffset + 2 ] );
1945
1946 T * tData = aData( p ).data() ;
1947
1948 index_t tDataOffset = 0 ;
1949
1950 int tCommTag = comm_tag( tMyRank, p ) + 1 ;
1951
1952 for ( index_t c : tChunkSizes )
1953 {
1954 comm_check( MPI_Irecv( &tData[ tDataOffset ],
1955 c,
1956 tCommType,
1957 p,
1958 tCommTag,
1959 gComm.world(),
1960 &tRequest[ tCount++ ] ) );
1961
1962 tDataOffset+= c;
1963 }
1964
1965 tOffset += 3 ;
1966 }
1967
1968 comm_check( MPI_Waitall( tCount, tRequest, tStatus ) );
1969
1970 free( tStatus );
1971 free( tRequest );
1972 free( tSizes );
1973
1974#endif
1975 }
1976
1977//==============================================================================
1978// STRINGS
1979//==============================================================================
1980
1982 void
1983 broadcast( Cell< string > & aData, const proc_t aRoot=0 );
1984
1986 void
1987 send( const string & aMessage, const proc_t aTarget=0 );
1988
1990 void
1991 receive( string & aMessage, const proc_t aSource=0 );
1992
1993//------------------------------------------------------------------------------
1994}
1995#endif //COMMTOOLS_HPP
BELFEM_ERROR
#define BELFEM_ERROR(aCheck,...)
Definition assert.hpp:264
BELFEM_ASSERT
#define BELFEM_ASSERT(aCheck,...)
Definition assert.hpp:244
cl_Communicator.hpp
gComm
belfem::Communicator gComm
Definition belfem.cpp:35
cl_Matrix.hpp
cl_Vector.hpp
belfem::Cell
A bounds-checked wrapper around std::vector.
Definition cl_Cell.hpp:45
belfem::Cell::set_size
void set_size(const size_t aSize)
New elements are value-initialized ( zero for arithmetic T, nullptr for pointers ).
Definition cl_Cell.hpp:199
belfem::Cell::size
size_t size() const
Definition cl_Cell.hpp:189
belfem::Cell::data
T * data()
Returns a borrowed pointer into contiguous storage, possibly nullptr when the Cell is empty ( test si...
Definition cl_Cell.hpp:129
belfem::Communicator
Global MPI communicator manager.
Definition cl_Communicator.hpp:60
belfem::Communicator::world
COMM_TYPE & world()
Definition cl_Communicator.cpp:251
belfem::Communicator::rank
proc_t rank() const
Definition cl_Communicator.hpp:192
belfem::Communicator::size
proc_t size() const
Definition cl_Communicator.hpp:200
belfem::Matrix
Dense column-major matrix.
Definition cl_BZ_Matrix.hpp:28
belfem::Matrix::spacing
size_t spacing() const
inter-column stride of the data container ( the leading dimension in BLAS terms; here always n_rows )
Definition cl_AR_Matrix.hpp:214
belfem::Matrix::capacity
size_t capacity() const
The logical element count n_rows() * n_cols() ( Armadillo's n_elem ).
Definition cl_AR_Matrix.hpp:202
belfem::Matrix::set_size
void set_size(const size_t aNumRows, const size_t aNumCols)
Definition cl_AR_Matrix.hpp:164
belfem::Matrix::n_rows
size_t n_rows() const
Definition cl_AR_Matrix.hpp:183
belfem::Matrix::n_cols
size_t n_cols() const
Definition cl_AR_Matrix.hpp:191
belfem::Matrix::data
T * data()
Returns a borrowed raw buffer under the layout contract stated in cl_Matrix.hpp; may be nullptr when ...
Definition cl_AR_Matrix.hpp:116
belfem::Vector
Column vector.
Definition cl_BZ_Vector.hpp:41
belfem::Vector::data
T * data()
Returns a borrowed raw buffer under the layout contract stated in cl_Vector.hpp; may be nullptr when ...
Definition cl_AR_Vector.hpp:163
belfem::Vector::set_size
void set_size(const size_t aNumRows)
change the size of the vector
Definition cl_AR_Vector.hpp:215
belfem::Vector::length
size_t length() const
get the length of the vector
Definition cl_AR_Vector.hpp:235
belfem::constant::c
const real c
speed of light in m/s ( exact ) http://physics.nist.gov/cgi-bin/cuu/Value?c
Definition constants.hpp:65
belfem
USER GUIDES:
Definition cl_Capacitor.cpp:16
belfem::comm_rank
proc_t comm_rank()
Definition commtools.cpp:30
belfem::comm_type
comm_t comm_type()
returns the MPI datatype handle for T
Definition commtypes.hpp:37
belfem::share
void share(Vector< T > &aData)
Sends only; it is not collective.
Definition commtools.hpp:1344
belfem::collect
void collect(Cell< T > &aData, const T aMyValue=0)
Receives only, one slot per rank.
Definition commtools.hpp:702
belfem::MaterialDependency::T
@ T
Definition cl_Material.hpp:125
belfem::receive
void receive(string &aMessage, const proc_t aSource)
Same contract as receive( T& ).
Definition commtools.cpp:348
belfem::comm_barrier
void comm_barrier()
Collective: every rank must reach it before any rank passes it.
Definition commtools.cpp:57
belfem::comm_t
int comm_t
Definition commtypes.hpp:27
belfem::comm_drain_check
void comm_drain_check(const char *aLabel)
Debug-build tripwire for the two-tags-per-pair ordering contract of the point-to-point fabric ( see c...
Definition commtools.cpp:117
belfem::broadcast
void broadcast(Cell< string > &aData, const proc_t aRoot)
Collective: every rank calls it, and aRoot's data lands on all ranks.
Definition commtools.cpp:276
belfem::proc_t
int proc_t
Definition commtypes.hpp:29
belfem::comm_splitcount
index_t comm_splitcount(const Vector< index_t > &aLengths, const proc_t aRoot)
Local, it does not communicate despite aRoot: the number of chunk messages comm_split() would produce...
Definition commtools.cpp:228
belfem::distribute
void distribute(Cell< T > &aData)
Sends only, not collective: aData has one entry per rank ( asserted ) and entry p goes to rank p; its...
Definition commtools.hpp:657
belfem::allreduce
void allreduce(const T *aSend, T *aRecv, const int aCount)
Collective MAX-reduction visible on every rank.
Definition commtools.hpp:203
belfem::comm_tag
int comm_tag(const proc_t aSource, const proc_t aTarget)
The base tag of an unordered rank pair.
Definition commtools.cpp:102
belfem::gMaxCommChunkLength
constexpr int gMaxCommChunkLength
Definition commtools.hpp:31
belfem::comm_split
Cell< int > comm_split(const index_t aLength)
Local, it does not communicate: splits aLength into chunks of at most gMaxCommChunkLength.
Definition commtools.cpp:207
belfem::index_t
uint32_t index_t
Definition typedefs.hpp:52
belfem::allreduce_min
void allreduce_min(const T *aSend, T *aRecv, const int aCount)
Collective MIN-reduction visible on every rank.
Definition commtools.hpp:245
belfem::comm_check
void comm_check(const int aErrorCode)
Any code other than MPI_SUCCESS aborts through BELFEM_ERROR ( release too ), which ends the whole run...
Definition commtools.cpp:38
belfem::send
void send(const string &aMessage, const proc_t aTarget)
Same contract as send( T ).
Definition commtools.cpp:340
belfem::comm_size
proc_t comm_size()
Definition commtools.cpp:22
corc.tests.test_curve_frame.p
p
Definition test_curve_frame.py:73
belfem::is_scalar
Definition commtools.hpp:38