test_mapweights_bcast.cpp File Reference

#include <iostream>
#include <vector>
#include "moab/Core.hpp"
#include "OfflineMap.h"
#include "moab/Remapping/TempestRemapper.hpp"
#include "moab/Remapping/TempestOnlineMap.hpp"
#include "moab/ParallelComm.hpp"
#include "Internals.hpp"
#include "TestRunner.hpp"
#include <netcdf.h>
#include <netcdf_par.h>

Include dependency graph for test_mapweights_bcast.cpp:

Go to the source code of this file.

Macros
#define	IS_BUILDING_MB
#define	NCERROREXIT(err, MSG)

Functions
void	test_tempest_map_bcast ()
void	read_buffered_map ()
void	read_map_from_disk ()
int	main (int argc, char **argv)

Variables
std::string	inMapFilename = ""

Macro Definition Documentation

IS_BUILDING_MB

#define IS_BUILDING_MB

Definition at line 24 of file test_mapweights_bcast.cpp.

NCERROREXIT

#define NCERROREXIT	(		err,
			MSG
	)

Value:

 if( err ) \
 { \
 std::cout << "Error (" << ( err ) << "): " << ( MSG ); \
 std::cout << "\nAborting with message: " << nc_strerror( err ) << "... \n"; \
 return; \
 }

Definition at line 62 of file test_mapweights_bcast.cpp.

Function Documentation

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 41 of file test_mapweights_bcast.cpp.

{
 if( argc > 1 )
 inMapFilename = std::string( argv[1] );
 else
 inMapFilename = TestDir + "unittest/outCS5ICOD5_map.nc";
 
 REGISTER_TEST( test_tempest_map_bcast );
 REGISTER_TEST( read_buffered_map );
 REGISTER_TEST( read_map_from_disk );
 
#ifdef MOAB_HAVE_MPI
 MPI_Init( &argc, &argv );
#endif
 int result = RUN_TESTS( argc, argv );
#ifdef MOAB_HAVE_MPI
 MPI_Finalize();
#endif
 return result;
}

References inMapFilename, read_buffered_map(), read_map_from_disk(), REGISTER_TEST, RUN_TESTS, and test_tempest_map_bcast().

read_buffered_map()

void read_buffered_map

(

)

Perform a series of checks to ensure that the local maps in each process still preserve map properties Let us run our usual tests to verify that the map data has been distributed correctly 1) Consistency check – row sum = 1.0 2) Disable conservation checks that require global column sums (and hence more communication/reduction) 3) Perform monotonicity check and ensure failures globally are same as serial case

Definition at line 70 of file test_mapweights_bcast.cpp.

{
 NcError error( NcError::silent_nonfatal );
 
 // Define our total buffer size to use
 NcFile* ncMap = NULL;
 NcVar *varRow = NULL, *varCol = NULL, *varS = NULL;
 
#ifdef MOAB_HAVE_MPI
 int mpierr = 0;
 const int rootProc = 0;
#endif
 // Main inputs;
 // 1. Map Filename
 // 2. Ownership info from the processes to be gathered on root
 const int nBufferSize = 64 * 1024; // 64 KB
 const int nNNZBytes = 2 * sizeof( int ) + sizeof( double );
 const int nMaxEntries = nBufferSize / nNNZBytes;
 int nA = 0, nB = 0, nVA = 0, nVB = 0, nS = 0;
 double dTolerance = 1e-08;
 
 int nprocs = 1, rank = 0;
#ifdef MOAB_HAVE_MPI
 MPI_Comm commW = MPI_COMM_WORLD;
 MPI_Comm_size( commW, &nprocs );
 MPI_Comm_rank( commW, &rank );
#endif
 
#ifdef MOAB_HAVE_MPI
 if( rank == rootProc )
#endif
 {
 ncMap = new NcFile( inMapFilename.c_str(), NcFile::ReadOnly );
 if( !ncMap->is_valid() )
 {
 _EXCEPTION1( "Unable to open input map file \"%s\"", inMapFilename.c_str() );
 }
 
 // Source and Target mesh resolutions
 NcDim* dimNA = ncMap->get_dim( "n_a" );
 if( dimNA == NULL )
 {
 _EXCEPTIONT( "Input map missing dimension \"n_a\"" );
 }
 else
 nA = dimNA->size();
 
 NcDim* dimNB = ncMap->get_dim( "n_b" );
 if( dimNB == NULL )
 {
 _EXCEPTIONT( "Input map missing dimension \"n_b\"" );
 }
 else
 nB = dimNB->size();
 
 NcDim* dimNVA = ncMap->get_dim( "nv_a" );
 if( dimNVA == NULL )
 {
 _EXCEPTIONT( "Input map missing dimension \"nv_a\"" );
 }
 else
 nVA = dimNVA->size();
 
 NcDim* dimNVB = ncMap->get_dim( "nv_b" );
 if( dimNVB == NULL )
 {
 _EXCEPTIONT( "Input map missing dimension \"nv_b\"" );
 }
 else
 nVB = dimNVB->size();
 
 // Read SparseMatrix entries
 NcDim* dimNS = ncMap->get_dim( "n_s" );
 if( dimNS == NULL )
 {
 _EXCEPTION1( "Map file \"%s\" does not contain dimension \"n_s\"", inMapFilename.c_str() );
 }
 else
 nS = dimNS->size(); // store total number of nonzeros
 
 varRow = ncMap->get_var( "row" );
 if( varRow == NULL )
 {
 _EXCEPTION1( "Map file \"%s\" does not contain variable \"row\"", inMapFilename.c_str() );
 }
 
 varCol = ncMap->get_var( "col" );
 if( varCol == NULL )
 {
 _EXCEPTION1( "Map file \"%s\" does not contain variable \"col\"", inMapFilename.c_str() );
 }
 
 varS = ncMap->get_var( "S" );
 if( varS == NULL )
 {
 _EXCEPTION1( "Map file \"%s\" does not contain variable \"S\"", inMapFilename.c_str() );
 }
 }
 
 // const int nTotalBytes = nS * nNNZBytes;
 int nEntriesRemaining = nS;
 int nBufferedReads = static_cast< int >( ceil( ( 1.0 * nS ) / ( 1.0 * nMaxEntries ) ) );
 int runData[6] = { nA, nB, nVA, nVB, nS, nBufferedReads };
 
 mpierr = MPI_Bcast( runData, 6, MPI_INT, rootProc, commW );
 CHECK_EQUAL( mpierr, 0 );
 
#ifdef MOAB_HAVE_MPI
 if( rank != rootProc )
 {
 nA = runData[0];
 nB = runData[1];
 nVA = runData[2];
 nVB = runData[3];
 nS = runData[4];
 nBufferedReads = runData[5];
 }
 else
 printf( "Global parameters: nA = %d, nB = %d, nS = %d\n", nA, nB, nS );
#else
 printf( "Global parameters: nA = %d, nB = %d, nS = %d\n", nA, nB, nS );
#endif
 
 std::vector< int > rowOwnership( nprocs );
 int nGRowPerPart = nB / nprocs;
 int nGRowRemainder = nB % nprocs; // Keep the remainder in root
 rowOwnership[0] = nGRowPerPart + nGRowRemainder;
 for( int ip = 1, roffset = rowOwnership[0]; ip < nprocs; ++ip )
 {
 roffset += nGRowPerPart;
 rowOwnership[ip] = roffset;
 }
 
 // Let us declare the map object for every process
 OfflineMap mapRemap;
 SparseMatrix< double >& sparseMatrix = mapRemap.GetSparseMatrix();
 
#ifdef MOAB_HAVE_MPI
 if( rank == rootProc )
#endif
 printf( "Parameters: nA=%d, nB=%d, nVA=%d, nVB=%d, nS=%d, nNNZBytes = %d, nBufferedReads = %d\n", nA, nB, nVA,
 nVB, nS, nNNZBytes, nBufferedReads );
 
 std::map< int, int > rowMap, colMap;
 int rindexMax = 0, cindexMax = 0;
 long offset = 0;
 
 /* Split the rows and send to processes in chunks
 Let us start the buffered read */
 for( int iRead = 0; iRead < nBufferedReads; ++iRead )
 {
 // pointers to the data
 DataArray1D< int > vecRow;
 DataArray1D< int > vecCol;
 DataArray1D< double > vecS;
 std::vector< std::vector< int > > dataPerProcess( nprocs );
 std::vector< int > nDataPerProcess( nprocs, 0 );
 for( int ip = 0; ip < nprocs; ++ip )
 dataPerProcess[ip].reserve( std::max( 100, nMaxEntries / nprocs ) );
 
#ifdef MOAB_HAVE_MPI
 if( rank == rootProc )
#endif
 {
 int nLocSize = std::min( nEntriesRemaining, static_cast< int >( ceil( nBufferSize * 1.0 / nNNZBytes ) ) );
 
 printf( "Reading file: elements %ld to %ld\n", offset, offset + nLocSize );
 
 // Allocate and resize based on local buffer size
 vecRow.Allocate( nLocSize );
 vecCol.Allocate( nLocSize );
 vecS.Allocate( nLocSize );
 
 varRow->set_cur( (long)( offset ) );
 varRow->get( &( vecRow[0] ), nLocSize );
 
 varCol->set_cur( (long)( offset ) );
 varCol->get( &( vecCol[0] ), nLocSize );
 
 varS->set_cur( (long)( offset ) );
 varS->get( &( vecS[0] ), nLocSize );
 
 // Decrement vecRow and vecCol
 for( size_t i = 0; i < vecRow.GetRows(); i++ )
 {
 vecRow[i]--;
 vecCol[i]--;
 
 // TODO: Fix this linear search to compute process ownership
 int pOwner = -1;
 if( rowOwnership[0] > vecRow[i] )
 pOwner = 0;
 else
 {
 for( int ip = 1; ip < nprocs; ++ip )
 {
 if( rowOwnership[ip - 1] <= vecRow[i] && rowOwnership[ip] > vecRow[i] )
 {
 pOwner = ip;
 break;
 }
 }
 }
 
 assert( pOwner >= 0 && pOwner < nprocs );
 dataPerProcess[pOwner].push_back( i );
 }
 
 offset += nLocSize;
 nEntriesRemaining -= nLocSize;
 
 for( int ip = 0; ip < nprocs; ++ip )
 nDataPerProcess[ip] = dataPerProcess[ip].size();
 }
 
 // Communicate the number of DoF data to expect from root
 int nEntriesComm = 0;
 mpierr = MPI_Scatter( nDataPerProcess.data(), 1, MPI_INT, &nEntriesComm, 1, MPI_INT, rootProc, commW );
 CHECK_EQUAL( mpierr, 0 );
 
#ifdef MOAB_HAVE_MPI
 if( rank == rootProc )
#endif
 {
 std::vector< MPI_Request > cRequests;
 std::vector< MPI_Status > cStats( 2 * nprocs - 2 );
 cRequests.reserve( 2 * nprocs - 2 );
 
 // First send out all the relevant data buffers to other process
 std::vector< std::vector< int > > dataRowColsAll( nprocs );
 std::vector< std::vector< double > > dataEntriesAll( nprocs );
 for( int ip = 1; ip < nprocs; ++ip )
 {
 const int nDPP = nDataPerProcess[ip];
 if( nDPP > 0 )
 {
 std::vector< int >& dataRowCols = dataRowColsAll[ip];
 dataRowCols.resize( 2 * nDPP );
 std::vector< double >& dataEntries = dataEntriesAll[ip];
 dataEntries.resize( nDPP );
 
 for( int ij = 0; ij < nDPP; ++ij )
 {
 dataRowCols[ij * 2] = vecRow[dataPerProcess[ip][ij]];
 dataRowCols[ij * 2 + 1] = vecCol[dataPerProcess[ip][ij]];
 dataEntries[ij] = vecS[dataPerProcess[ip][ij]];
 }
 
 MPI_Request rcsend, dsend;
 mpierr = MPI_Isend( dataRowCols.data(), 2 * nDPP, MPI_INT, ip, iRead * 1000, commW, &rcsend );
 CHECK_EQUAL( mpierr, 0 );
 mpierr = MPI_Isend( dataEntries.data(), nDPP, MPI_DOUBLE, ip, iRead * 1000 + 1, commW, &dsend );
 CHECK_EQUAL( mpierr, 0 );
 
 cRequests.push_back( rcsend );
 cRequests.push_back( dsend );
 }
 }
 
 // Next perform all necessary local work while we wait for the buffers to be sent out
 // Compute an offset for the rows and columns by creating a local to global mapping for rootProc
 int rindex = 0, cindex = 0;
 assert( dataPerProcess[0].size() - nEntriesComm == 0 ); // sanity check
 for( int i = 0; i < nEntriesComm; ++i )
 {
 const int& vecRowValue = vecRow[dataPerProcess[0][i]];
 const int& vecColValue = vecCol[dataPerProcess[0][i]];
 
 std::map< int, int >::iterator riter = rowMap.find( vecRowValue );
 if( riter == rowMap.end() )
 {
 rowMap[vecRowValue] = rindexMax;
 rindex = rindexMax;
 rindexMax++;
 }
 else
 rindex = riter->second;
 
 std::map< int, int >::iterator citer = colMap.find( vecColValue );
 if( citer == colMap.end() )
 {
 colMap[vecColValue] = cindexMax;
 cindex = cindexMax;
 cindexMax++;
 }
 else
 cindex = citer->second;
 
 sparseMatrix( rindex, cindex ) = vecS[i];
 }
 
 // Wait until all communication is pushed out
 mpierr = MPI_Waitall( cRequests.size(), cRequests.data(), cStats.data() );
 CHECK_EQUAL( mpierr, 0 );
 } // if( rank == rootProc )
#ifdef MOAB_HAVE_MPI
 else
 {
 if( nEntriesComm > 0 )
 {
 MPI_Request cRequests[2];
 MPI_Status cStats[2];
 /* code */
 // create buffers to receive the data from root process
 std::vector< int > dataRowCols( 2 * nEntriesComm );
 vecRow.Allocate( nEntriesComm );
 vecCol.Allocate( nEntriesComm );
 vecS.Allocate( nEntriesComm );
 
 /* TODO: combine row and column scatters together. Save one communication by better packing */
 // Scatter the rows, cols and entries to the processes according to the partition
 mpierr = MPI_Irecv( dataRowCols.data(), 2 * nEntriesComm, MPI_INT, rootProc, iRead * 1000, commW,
 &cRequests[0] );
 CHECK_EQUAL( mpierr, 0 );
 
 mpierr = MPI_Irecv( vecS, nEntriesComm, MPI_DOUBLE, rootProc, iRead * 1000 + 1, commW, &cRequests[1] );
 CHECK_EQUAL( mpierr, 0 );
 
 // Wait until all communication is pushed out
 mpierr = MPI_Waitall( 2, cRequests, cStats );
 CHECK_EQUAL( mpierr, 0 );
 
 // Compute an offset for the rows and columns by creating a local to global mapping
 int rindex = 0, cindex = 0;
 for( int i = 0; i < nEntriesComm; ++i )
 {
 const int& vecRowValue = dataRowCols[i * 2];
 const int& vecColValue = dataRowCols[i * 2 + 1];
 
 std::map< int, int >::iterator riter = rowMap.find( vecRowValue );
 if( riter == rowMap.end() )
 {
 rowMap[vecRowValue] = rindexMax;
 rindex = rindexMax;
 rindexMax++;
 }
 else
 rindex = riter->second;
 vecRow[i] = rindex;
 
 std::map< int, int >::iterator citer = colMap.find( vecColValue );
 if( citer == colMap.end() )
 {
 colMap[vecColValue] = cindexMax;
 cindex = cindexMax;
 cindexMax++;
 }
 else
 cindex = citer->second;
 vecCol[i] = cindex;
 
 sparseMatrix( rindex, cindex ) = vecS[i];
 }
 
 // clear the local buffer
 dataRowCols.clear();
 } // if( nEntriesComm > 0 )
 } // if( rank != rootProc )
 
 MPI_Barrier( commW );
#endif
 }
 
#ifdef MOAB_HAVE_MPI
 if( rank == rootProc )
#endif
 {
 assert( nEntriesRemaining == 0 );
 ncMap->close();
 }
 
 /**
 * Perform a series of checks to ensure that the local maps in each process still preserve map properties
 * Let us run our usual tests to verify that the map data has been distributed correctly
 * 1) Consistency check -- row sum = 1.0
 * 2) Disable conservation checks that require global column sums (and hence more communication/reduction)
 * 3) Perform monotonicity check and ensure failures globally are same as serial case
 */
 printf( "Rank %d: sparsematrix size = %d x %d\n", rank, sparseMatrix.GetRows(), sparseMatrix.GetColumns() );
 
 // consistency: row sums
 {
 int isConsistentP = mapRemap.IsConsistent( dTolerance );
 if( 0 != isConsistentP )
 {
 printf( "Rank %d failed consistency checks...\n", rank );
 }
 CHECK_EQUAL( 0, isConsistentP );
 }
 
 // conservation: we will disable this conservation check for now.
 // int isConservativeP = mapRemap.IsConservative( dTolerance );
 // CHECK_EQUAL( 0, isConservativeP );
 
 // monotonicity: local failures
 {
 int isMonotoneP = mapRemap.IsMonotone( dTolerance );
 // Accumulate sum of failures to ensure that it is exactly same as serial case
 int isMonotoneG;
 mpierr = MPI_Allreduce( &isMonotoneP, &isMonotoneG, 1, MPI_INT, MPI_SUM, commW );
 CHECK_EQUAL( mpierr, 0 );
 
 // 4600 monotonicity fails in serial and parallel for the test map file
 CHECK_EQUAL( 4600, isMonotoneG );
 }
 
 delete ncMap;
 
 return;
}

References CHECK_EQUAL, moab::error(), inMapFilename, MPI_COMM_WORLD, rank, and size.

Referenced by main().

read_map_from_disk()

void read_map_from_disk

(

)

Definition at line 481 of file test_mapweights_bcast.cpp.

{
 moab::ErrorCode rval;
 NcError error( NcError::verbose_nonfatal );
 MPI_Comm commW = MPI_COMM_WORLD;
 std::vector< int > tmp_owned_ids;
 double dTolerance = 1e-08;
 
 std::string remap_weights_filename = TestDir + "unittest/outCS5ICOD5_map.nc";
 
 moab::Interface* mb = new moab::Core();
 moab::ParallelComm pcomm( mb, commW );
 
 moab::TempestRemapper remapper( mb, &pcomm );
 remapper.meshValidate = true;
 remapper.constructEdgeMap = true;
 
 // Do not create new filesets; Use the sets from our respective applications
 remapper.initialize( false );
 
 moab::TempestOnlineMap onlinemap( &remapper );
 
 rval = onlinemap.ReadParallelMap( remap_weights_filename.c_str(), tmp_owned_ids, true );
 CHECK_EQUAL( rval, moab::MB_SUCCESS );
 
 // consistency: row sums
 int isConsistentP = onlinemap.IsConsistent( dTolerance );
 CHECK_EQUAL( 0, isConsistentP );
 
 // conservation: we will disable this conservation check for now
 // since we do not have information about source and target areas
 // to accurately decipher conservation data satisfaction
 // int isConservativeP = onlinemap.IsConservative( dTolerance );
 // CHECK_EQUAL( 0, isConservativeP );
 
 // monotonicity: global failures
 int isMonotoneP = onlinemap.IsMonotone( dTolerance );
 // Accumulate sum of failures to ensure that it is exactly same as serial case
 // 4600 fails in serial. What about parallel ? Check and confirm.
 CHECK_EQUAL( 4600, isMonotoneP );
}

References CHECK_EQUAL, moab::TempestRemapper::constructEdgeMap, moab::error(), ErrorCode, moab::TempestRemapper::initialize(), moab::TempestOnlineMap::IsConsistent(), moab::TempestOnlineMap::IsMonotone(), mb, MB_SUCCESS, moab::TempestRemapper::meshValidate, MPI_COMM_WORLD, and moab::TempestOnlineMap::ReadParallelMap().

Referenced by main().

test_tempest_map_bcast()

void test_tempest_map_bcast

(

)

Perform a series of checks to ensure that the local maps in each process still preserve map properties Let us run our usual tests to verify that the map data has been distributed correctly 1) Consistency check – row sum = 1.0 2) Disable conservation checks that require global column sums (and hence more communication/reduction) 3) Perform monotonicity check and ensure failures globally are same as serial case

Definition at line 523 of file test_mapweights_bcast.cpp.

{
 /*
 *
 * Load mapping file generated with CS5 and ICOD5 meshes with FV-FV
 * projection for field data. It has following attributes:
 *
 * Analyzing map
 * ..Per-dof consistency (tol 1.00000e-08).. PASS
 * ..Per-dof conservation (tol 1.00000e-08).. PASS
 * ..Weights within range [-10,+10].. Done
 * ..
 * .. Total nonzero entries: 8976
 * .. Column index min/max: 1 / 150 (150 source dofs)
 * .. Row index min/max: 1 / 252 (252 target dofs)
 * .. Source area / 4pi: 9.999999999991666e-01
 * .. Source area min/max: 7.758193626656797e-02 / 9.789674024202324e-02
 * .. Target area / 4pi: 9.999999999999928e-01
 * .. Target area min/max: 3.418848585325412e-02 / 5.362041648788460e-02
 * .. Source frac min/max: 9.999999999997851e-01 / 1.000000000003865e+00
 * .. Target frac min/max: 9.999999999993099e-01 / 1.000000000000784e+00
 * .. Map weights min/max: -2.062659472710135e-01 / 1.143815352351317e+00
 * .. Row sums min/max: 9.999999999993099e-01 / 1.000000000000784e+00
 * .. Consist. err min/max: -6.901146321069973e-13 / 7.838174553853605e-13
 * .. Col wt.sums min/max: 9.999999999997854e-01 / 1.000000000003865e+00
 * .. Conserv. err min/max: -2.146061106600428e-13 / 3.865352482534945e-12
 * ..
 * ..Histogram of nonzero entries in sparse matrix
 * ....Column 1: Number of nonzero entries (bin minimum)
 * ....Column 2: Number of columns with that many nonzero values
 * ....Column 3: Number of rows with that many nonzero values
 * ..[[25, 0, 2],[29, 0, 8],[30, 0, 10],[31, 150, 232]]
 * ..
 * ..Histogram of weights
 * ....Column 1: Lower bound on weights
 * ....Column 2: Upper bound on weights
 * ....Column 3: # of weights in that bin
 * ..[[-1, 0, 4577],[0, 1, 4365],[1, 2, 34]]
 *
 */
 NcError error( NcError::silent_nonfatal );
 MPI_Comm commW = MPI_COMM_WORLD;
 const int rootProc = 0;
 const std::string outFilename = inMapFilename;
 double dTolerance = 1e-08;
 int mpierr = 0;
 
 int nprocs, rank;
 MPI_Comm_size( MPI_COMM_WORLD, &nprocs );
 MPI_Comm_rank( MPI_COMM_WORLD, &rank );
 
 const double dFillValueOverride = 0.0;
 OfflineMap mapRemap;
 if( rank == 0 )
 {
 std::cout << "Reading file: " << outFilename << std::endl;
 mapRemap.Read( outFilename );
 mapRemap.SetFillValueOverrideDbl( dFillValueOverride );
 mapRemap.SetFillValueOverride( static_cast< float >( dFillValueOverride ) );
 
 int isConsistent = mapRemap.IsConsistent( dTolerance );
 CHECK_EQUAL( 0, isConsistent );
 int isConservative = mapRemap.IsConservative( dTolerance );
 CHECK_EQUAL( 0, isConservative );
 int isMonotone = mapRemap.IsMonotone( dTolerance );
 CHECK_EQUAL( 4600, isMonotone );
 
 // verify that the source and target mesh areas are equal
 DataArray1D< double >& srcAreas = mapRemap.GetSourceAreas();
 DataArray1D< double >& tgtAreas = mapRemap.GetTargetAreas();
 
 double totSrcArea = 0.0;
 for( size_t i = 0; i < srcAreas.GetRows(); ++i )
 totSrcArea += srcAreas[i];
 
 double totTgtArea = 0.0;
 for( size_t i = 0; i < tgtAreas.GetRows(); ++i )
 totTgtArea += tgtAreas[i];
 
 CHECK_REAL_EQUAL( totSrcArea, totTgtArea, 1e-10 );
 }
 
 MPI_Barrier( commW );
 
 SparseMatrix< double >& sparseMatrix = mapRemap.GetSparseMatrix();
 
 // Get the row, col and entry data tupl
 // Retrieve all data from the map operator on root process
 DataArray1D< int > dataRowsGlobal, dataColsGlobal;
 DataArray1D< double > dataEntriesGlobal;
 if( rank == 0 )
 {
 sparseMatrix.GetEntries( dataRowsGlobal, dataColsGlobal, dataEntriesGlobal );
 std::cout << "Total NNZ in remap matrix = " << dataRowsGlobal.GetRows() << std::endl;
 }
 
 int *rg = nullptr, *cg = nullptr;
 double* de = nullptr;
 int nMapGlobalRowCols[3] = { sparseMatrix.GetRows(), sparseMatrix.GetColumns(),
 static_cast< int >( dataEntriesGlobal.GetRows() ) };
 
 /* split the rows and send to processes in chunks */
 std::vector< int > rowAllocation;
 std::vector< int > sendCounts;
 std::vector< int > displs;
 const int NDATA = 3;
 if( rank == 0 )
 {
 rg = dataRowsGlobal;
 cg = dataColsGlobal;
 de = dataEntriesGlobal;
 
 // Let us do a trivial partitioning of the row data
 rowAllocation.resize( nprocs * NDATA );
 displs.resize( nprocs );
 sendCounts.resize( nprocs );
 
 int nRowPerPart = nMapGlobalRowCols[0] / nprocs;
 int remainder = nMapGlobalRowCols[0] % nprocs; // Keep the remainder in root
 // printf( "nRowPerPart = %d, remainder = %d\n", nRowPerPart, remainder );
 int sum = 0;
 for( int i = 0; i < nprocs; ++i )
 {
 rowAllocation[i * NDATA] = nRowPerPart + ( i == 0 ? remainder : 0 );
 rowAllocation[i * NDATA + 1] = std::find( rg, rg + dataRowsGlobal.GetRows(), sum ) - rg;
 sum += rowAllocation[i * NDATA];
 displs[i] = rowAllocation[i * NDATA + 1];
 }
 for( int i = 0; i < nprocs - 1; ++i )
 {
 rowAllocation[i * NDATA + 2] = rowAllocation[( i + 1 ) * NDATA + 1] - rowAllocation[i * NDATA + 1];
 sendCounts[i] = rowAllocation[i * NDATA + 2];
 }
 rowAllocation[( nprocs - 1 ) * NDATA + 2] = nMapGlobalRowCols[2] - displs[nprocs - 1];
 sendCounts[nprocs - 1] = rowAllocation[( nprocs - 1 ) * NDATA + 2];
 
 // for( int i = 0; i < nprocs; i++ )
 // printf( "sendcounts[%d] = %d, %d, %d\tdispls[%d] = %d, %d\n", i, rowAllocation[i * NDATA],
 // rowAllocation[i * NDATA + 1], rowAllocation[i * NDATA + 2], i, displs[i], sendCounts[i] );
 }
 
 mpierr = MPI_Bcast( nMapGlobalRowCols, 3, MPI_INT, rootProc, commW );
 CHECK_EQUAL( mpierr, 0 );
 
 int allocationSize[NDATA] = { 0, 0, 0 };
 mpierr = MPI_Scatter( rowAllocation.data(), NDATA, MPI_INT, &allocationSize, NDATA, MPI_INT, rootProc, commW );
 CHECK_EQUAL( mpierr, 0 );
 
 // create buffers to receive the data from root process
 DataArray1D< int > dataRows, dataCols;
 DataArray1D< double > dataEntries;
 dataRows.Allocate( allocationSize[2] );
 dataCols.Allocate( allocationSize[2] );
 dataEntries.Allocate( allocationSize[2] );
 
 /* TODO: combine row and column scatters together. Save one communication by better packing */
 // Scatter the rows, cols and entries to the processes according to the partition
 mpierr = MPI_Scatterv( rg, sendCounts.data(), displs.data(), MPI_INT, (int*)( dataRows ), dataRows.GetByteSize(),
 MPI_INT, rootProc, commW );
 CHECK_EQUAL( mpierr, 0 );
 
 mpierr = MPI_Scatterv( cg, sendCounts.data(), displs.data(), MPI_INT, (int*)( dataCols ), dataCols.GetByteSize(),
 MPI_INT, rootProc, commW );
 CHECK_EQUAL( mpierr, 0 );
 
 mpierr = MPI_Scatterv( de, sendCounts.data(), displs.data(), MPI_DOUBLE, (double*)( dataEntries ),
 dataEntries.GetByteSize(), MPI_DOUBLE, rootProc, commW );
 CHECK_EQUAL( mpierr, 0 );
 
 // Compute an offset for the rows and columns by creating a local to global mapping
 std::vector< int > rowMap, colMap;
 for( int i = 0; i < allocationSize[2]; ++i )
 {
 int rindex, cindex;
 std::vector< int >::iterator riter = std::find( rowMap.begin(), rowMap.end(), dataRows[i] );
 if( riter == rowMap.end() )
 {
 rowMap.push_back( dataRows[i] );
 rindex = rowMap.size() - 1;
 }
 else
 rindex = riter - rowMap.begin();
 dataRows[i] = rindex;
 
 std::vector< int >::iterator citer = std::find( colMap.begin(), colMap.end(), dataCols[i] );
 if( citer == colMap.end() )
 {
 colMap.push_back( dataCols[i] );
 cindex = colMap.size() - 1;
 }
 else
 cindex = citer - colMap.begin();
 dataCols[i] = cindex;
 }
 
 // Now set the data received from root onto the sparse matrix
 sparseMatrix.SetEntries( dataRows, dataCols, dataEntries );
 
 /**
 * Perform a series of checks to ensure that the local maps in each process still preserve map properties
 * Let us run our usual tests to verify that the map data has been distributed correctly
 * 1) Consistency check -- row sum = 1.0
 * 2) Disable conservation checks that require global column sums (and hence more communication/reduction)
 * 3) Perform monotonicity check and ensure failures globally are same as serial case
 */
 
 // consistency: row sums
 int isConsistentP = mapRemap.IsConsistent( dTolerance );
 CHECK_EQUAL( 0, isConsistentP );
 
 // conservation: we will disable this conservation check for now.
 // int isConservativeP = mapRemap.IsConservative( dTolerance );
 // CHECK_EQUAL( 0, isConservativeP );
 
 // monotonicity: local failures
 int isMonotoneP = mapRemap.IsMonotone( dTolerance );
 // Accumulate sum of failures to ensure that it is exactly same as serial case
 int isMonotoneG;
 mpierr = MPI_Allreduce( &isMonotoneP, &isMonotoneG, 1, MPI_INT, MPI_SUM, commW );
 CHECK_EQUAL( mpierr, 0 );
 // 4600 fails in serial. What about parallel ? Check and confirm.
 CHECK_EQUAL( 4600, isMonotoneG );
 
 return;
}

References CHECK_EQUAL, CHECK_REAL_EQUAL, moab::error(), inMapFilename, MPI_COMM_WORLD, rank, and moab::sum().

Referenced by main().

Variable Documentation

inMapFilename

std::string inMapFilename = ""

Definition at line 35 of file test_mapweights_bcast.cpp.

Referenced by main(), read_buffered_map(), and test_tempest_map_bcast().