|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
#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.
Defines | |
| #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 = "" |
| #define IS_BUILDING_MB |
Definition at line 24 of file test_mapweights_bcast.cpp.
| #define NCERROREXIT | ( | err, | |
| MSG | |||
| ) |
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.
| int main | ( | int | argc, |
| char ** | argv | ||
| ) |
Definition at line 41 of file test_mapweights_bcast.cpp.
References inMapFilename, read_buffered_map(), read_map_from_disk(), REGISTER_TEST, PointInVol::result, RUN_TESTS, and test_tempest_map_bcast().
{
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;
}
| 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.
References CHECK_EQUAL, moab::error(), inMapFilename, MPI_COMM_WORLD, rank, and size.
Referenced by main().
{
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;
}
| void read_map_from_disk | ( | ) |
Definition at line 481 of file test_mapweights_bcast.cpp.
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().
{
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 );
}
| 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.
References CHECK_EQUAL, CHECK_REAL_EQUAL, moab::error(), inMapFilename, MPI_COMM_WORLD, rank, and moab::sum().
Referenced by main().
{
/*
*
* 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;
}
| std::string inMapFilename = "" |
Definition at line 35 of file test_mapweights_bcast.cpp.
Referenced by main(), read_buffered_map(), and test_tempest_map_bcast().
1.7.6.1