Skip to content

Commit

Permalink
Tpetra:
Browse files Browse the repository at this point in the history 
Specialization of sortCrsEntries and sortAndMergeCrsEntries that don't
sort/merge CRS values.  These procedures will be used by Tpetra::CrsGraph.

@trilinos/tpetra

Part of: trilinos#2267
  • Loading branch information
@tjfulle
tjfulle committed Mar 8, 2018
1 parent 541e09c commit d482822
Show file tree
Hide file tree
Showing 2 changed files with 193 additions and 15 deletions.
85 changes: 70 additions & 15 deletions packages/tpetra/core/src/Tpetra_Import_Util2.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,13 +75,22 @@ sortCrsEntries (const Teuchos::ArrayView<size_t>& CRS_rowptr,
const Teuchos::ArrayView<Ordinal>& CRS_colind,
const Teuchos::ArrayView<Scalar>&CRS_vals);

template<typename Ordinal>
void
sortCrsEntries (const Teuchos::ArrayView<size_t>& CRS_rowptr,
const Teuchos::ArrayView<Ordinal>& CRS_colind);

template<typename rowptr_array_type, typename colind_array_type, typename vals_array_type>
void
sortCrsEntries (const rowptr_array_type& CRS_rowptr,
const colind_array_type& CRS_colind,
const vals_array_type& CRS_vals);

template<typename rowptr_array_type, typename colind_array_type>
void
sortCrsEntries (const rowptr_array_type& CRS_rowptr,
const colind_array_type& CRS_colind);

/// \brief Sort and merge the entries of the (raw CSR) matrix by
/// column index within each row.
///
Expand All @@ -92,6 +101,11 @@ sortAndMergeCrsEntries (const Teuchos::ArrayView<size_t>& CRS_rowptr,
const Teuchos::ArrayView<Ordinal>& CRS_colind,
const Teuchos::ArrayView<Scalar>& CRS_vals);

template<typename Ordinal>
void
sortAndMergeCrsEntries (const Teuchos::ArrayView<size_t>& CRS_rowptr,
const Teuchos::ArrayView<Ordinal>& CRS_colind);

/// \brief lowCommunicationMakeColMapAndReindex
///
/// If you know the owning PIDs already, you can make the colmap a lot
Expand Down Expand Up @@ -168,9 +182,11 @@ sortCrsEntries (const Teuchos::ArrayView<size_t> &CRS_rowptr,
size_t start=CRS_rowptr[i];
if(start >= nnz) continue;

Scalar* locValues = &CRS_vals[start];
size_t NumEntries = CRS_rowptr[i+1] - start;
Ordinal* locIndices = &CRS_colind[start];
Teuchos::ArrayRCP<Scalar> locValues;
if (CRS_vals.size())
locValues = Teuchos::arcp<Scalar>(&CRS_vals[start], 0, NumEntries, false);
Teuchos::ArrayRCP<Ordinal> locIndices(&CRS_colind[start], 0, NumEntries, false);

Ordinal n = NumEntries;
Ordinal m = 1;
Expand All @@ -183,9 +199,11 @@ sortCrsEntries (const Teuchos::ArrayView<size_t> &CRS_rowptr,
for(Ordinal k = j; k >= 0; k-=m) {
if(locIndices[k+m] >= locIndices[k])
break;
Scalar dtemp = locValues[k+m];
locValues[k+m] = locValues[k];
locValues[k] = dtemp;
if (!Teuchos::is_null(locValues)) {
Scalar dtemp = locValues[k+m];
locValues[k+m] = locValues[k];
locValues[k] = dtemp;
}
Ordinal itemp = locIndices[k+m];
locIndices[k+m] = locIndices[k];
locIndices[k] = itemp;
Expand All @@ -196,6 +214,16 @@ sortCrsEntries (const Teuchos::ArrayView<size_t> &CRS_rowptr,
}
}

template<typename Ordinal>
void
sortCrsEntries (const Teuchos::ArrayView<size_t> &CRS_rowptr,
const Teuchos::ArrayView<Ordinal> & CRS_colind)
{
// Generate dummy values array
Teuchos::ArrayView<Tpetra::Details::DefaultTypes::scalar_type> CRS_vals;
sortCrsEntries (CRS_rowptr, CRS_colind, CRS_vals);
}

namespace Impl {

template<class RowOffsetsType, class ColumnIndicesType, class ValuesType>
Expand Down Expand Up @@ -238,9 +266,11 @@ class SortCrsEntries {
if (ind_(sk+m) >= ind_(sk)) {
break;
}
const scalar_type dtemp = val_(sk+m);
val_(sk+m) = val_(sk);
val_(sk) = dtemp;
if (val_.dimension_0()) {
const scalar_type dtemp = val_(sk+m);
val_(sk+m) = val_(sk);
val_(sk) = dtemp;
}
const ordinal_type itemp = ind_(sk+m);
ind_(sk+m) = ind_(sk);
ind_(sk) = itemp;
Expand Down Expand Up @@ -293,6 +323,17 @@ sortCrsEntries (const rowptr_array_type& CRS_rowptr,
vals_array_type>::sortCrsEntries (CRS_rowptr, CRS_colind, CRS_vals);
}

template<typename rowptr_array_type, typename colind_array_type>
void
sortCrsEntries (const rowptr_array_type& CRS_rowptr,
const colind_array_type& CRS_colind)
{
// Generate dummy values array
colind_array_type CRS_vals;
Impl::SortCrsEntries<rowptr_array_type, colind_array_type,
colind_array_type>::sortCrsEntries (CRS_rowptr, CRS_colind, CRS_vals);
}

// Note: This should get merged with the other Tpetra sort routines eventually.
template<typename Scalar, typename Ordinal>
void
Expand All @@ -318,9 +359,11 @@ sortAndMergeCrsEntries (const Teuchos::ArrayView<size_t> &CRS_rowptr,
CRS_rowptr[i] = old_curr;
if(old_rowptr_i >= nnz) continue;

Scalar* locValues = &CRS_vals[old_rowptr_i];
size_t NumEntries = CRS_rowptr[i+1] - old_rowptr_i;
Ordinal* locIndices = &CRS_colind[old_rowptr_i];
Teuchos::ArrayRCP<Scalar> locValues;
if (CRS_vals.size())
locValues = Teuchos::arcp<Scalar>(&CRS_vals[old_rowptr_i], 0, NumEntries, false);
Teuchos::ArrayRCP<Ordinal> locIndices(&CRS_colind[old_rowptr_i], 0, NumEntries, false);

// Sort phase
Ordinal n = NumEntries;
Expand All @@ -332,9 +375,11 @@ sortAndMergeCrsEntries (const Teuchos::ArrayView<size_t> &CRS_rowptr,
for(Ordinal k = j; k >= 0; k-=m) {
if(locIndices[k+m] >= locIndices[k])
break;
Scalar dtemp = locValues[k+m];
locValues[k+m] = locValues[k];
locValues[k] = dtemp;
if (!Teuchos::is_null(locValues)) {
Scalar dtemp = locValues[k+m];
locValues[k+m] = locValues[k];
locValues[k] = dtemp;
}
Ordinal itemp = locIndices[k+m];
locIndices[k+m] = locIndices[k];
locIndices[k] = itemp;
Expand All @@ -346,14 +391,14 @@ sortAndMergeCrsEntries (const Teuchos::ArrayView<size_t> &CRS_rowptr,
// Merge & shrink
for(size_t j=old_rowptr_i; j < CRS_rowptr[i+1]; j++) {
if(j > old_rowptr_i && CRS_colind[j]==CRS_colind[new_curr-1]) {
CRS_vals[new_curr-1] += CRS_vals[j];
if (CRS_vals.size()) CRS_vals[new_curr-1] += CRS_vals[j];
}
else if(new_curr==j) {
new_curr++;
}
else {
CRS_colind[new_curr] = CRS_colind[j];
CRS_vals[new_curr] = CRS_vals[j];
if (CRS_vals.size()) CRS_vals[new_curr] = CRS_vals[j];
new_curr++;
}
}
Expand All @@ -363,6 +408,16 @@ sortAndMergeCrsEntries (const Teuchos::ArrayView<size_t> &CRS_rowptr,
CRS_rowptr[NumRows] = new_curr;
}

template<typename Ordinal>
void
sortAndMergeCrsEntries (const Teuchos::ArrayView<size_t> &CRS_rowptr,
const Teuchos::ArrayView<Ordinal> & CRS_colind)
{
Teuchos::ArrayView<Tpetra::Details::DefaultTypes::scalar_type> CRS_vals;
return sortAndMergeCrsEntries(CRS_rowptr, CRS_colind, CRS_vals);
}


template <typename LocalOrdinal, typename GlobalOrdinal, typename Node>
void
lowCommunicationMakeColMapAndReindex (const Teuchos::ArrayView<const size_t> &rowptr,
Expand Down
123 changes: 123 additions & 0 deletions packages/tpetra/core/test/ImportExport2/ImportExport2_UnitTests.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,8 @@

#include <map>
#include <vector>
#include <random>
#include <algorithm>
#include <Teuchos_OrdinalTraits.hpp>
#include <Teuchos_ScalarTraits.hpp>
#include <Teuchos_VerboseObject.hpp>
Expand Down Expand Up @@ -1988,6 +1990,9 @@ TEUCHOS_UNIT_TEST_TEMPLATE_3_DECL( Import_Util, UnpackAndCombineWithOwningPIDs,
colind[i]=Bmap.getLocalElement(colind[i]);
}

// Make copies to test non-value interface to sortCrsEntries
Teuchos::Array<size_t> rowptr2(rowptr);

// Sort the GIDs
Tpetra::Import_Util::sortCrsEntries<Scalar, GO> (rowptr (), colind (), vals ());

Expand All @@ -2008,6 +2013,123 @@ TEUCHOS_UNIT_TEST_TEMPLATE_3_DECL( Import_Util, UnpackAndCombineWithOwningPIDs,
}


// Unit Test the functionality in Tpetra_Import_Util
TEUCHOS_UNIT_TEST_TEMPLATE_3_DECL( Import_Util, SortCrsEntries, LO, GO, Scalar) {

typedef typename Teuchos::Array<GO>::size_type size_type;

int max_num_entries_per_row = 7; // should be odd
int num_cols = 15; // should be odd

Teuchos::Array<size_t> rowptr, rowptr2;
Teuchos::Array<GO> colind, colind2;
Teuchos::Array<Scalar> vals, vals2;

// Fill the CRS arrays, use random values
std::random_device rd;
std::mt19937 gen(rd());
std::uniform_real_distribution<Scalar> dist(1.0, 2.0);
int row = 0;
while (true) {
int m = (max_num_entries_per_row - 1) / 2;
int start = std::max(0, row - m);
int num_cols_this_row = std::min(row + 1 + m,
std::min(max_num_entries_per_row, m+num_cols-row));
int end = start + num_cols_this_row;

rowptr.push_back(static_cast<size_t>(colind.size()));
rowptr2.push_back(static_cast<size_t>(colind2.size()));

for (int col=start; col<end; col++) {
colind.push_back(col);
colind2.push_back(col);

Scalar rands = dist(gen);
vals.push_back(rands);
vals2.push_back(rands);

// Create duplicate for merge test
colind2.push_back(col);
vals2.push_back(rands);
}

if (row > 1 && num_cols_this_row == 1 + m) break;
row++;

}
size_type num_entries = colind.size();
rowptr.push_back(static_cast<size_t>(num_entries));

size_type num_entries2 = colind2.size();
rowptr2.push_back(static_cast<size_t>(num_entries2));

// Randomly shuffle values and column indices
Teuchos::Array<Scalar> vals_rand(num_entries);
Teuchos::Array<GO> colind_rand1(num_entries), colind_rand2(num_entries);
for (size_t i=0; i < static_cast<size_t>(rowptr.size()-1); ++i) {
size_t num_cols_this_row = rowptr[i+1] - rowptr[i];
Teuchos::Array<GO> ix(num_cols_this_row);
std::iota(ix.begin(), ix.end(), rowptr[i]);
std::shuffle(ix.begin(), ix.end(), gen);
for (size_t j=rowptr[i], k=0; j < rowptr[i+1]; ++j, ++k) {
vals_rand[j] = vals[ix[k]];
colind_rand1[j] = colind[ix[k]];
colind_rand2[j] = colind_rand1[j];
}
}

// Sort the GIDs and associated values
Tpetra::Import_Util::sortCrsEntries<Scalar, GO>(rowptr(), colind_rand1(), vals_rand());

// Sort the GIDs
Tpetra::Import_Util::sortCrsEntries<GO>(rowptr(), colind_rand2());

// Compare
TEST_COMPARE_ARRAYS(colind, colind_rand1);
TEST_COMPARE_ARRAYS(colind, colind_rand2);
TEST_COMPARE_FLOATING_ARRAYS(vals, vals_rand, 1.e-12);

// Randomly shuffle values and column indices for merge test
Teuchos::Array<Scalar> vals2_rand(num_entries2);
Teuchos::Array<GO> colind2_rand1(num_entries2), colind2_rand2(num_entries2);
for (size_t i=0; i < static_cast<size_t>(rowptr2.size()-1); ++i) {
size_t num_cols_this_row = rowptr2[i+1] - rowptr2[i];
Teuchos::Array<GO> ix(num_cols_this_row);
std::iota(ix.begin(), ix.end(), rowptr2[i]);
std::shuffle(ix.begin(), ix.end(), gen);
for (size_t j=rowptr2[i], k=0; j < rowptr2[i+1]; ++j, ++k) {
vals2_rand[j] = vals2[ix[k]];
colind2_rand1[j] = colind2[ix[k]];
colind2_rand2[j] = colind2_rand1[j];
}
}

// Sort and merge the GIDs and associated values
Teuchos::Array<size_t> rowptr2_1(rowptr2);
Tpetra::Import_Util::sortAndMergeCrsEntries<Scalar, GO>(rowptr2(), colind2_rand1(), vals2_rand());

// Sort and merge the GIDs
Tpetra::Import_Util::sortAndMergeCrsEntries<GO>(rowptr2_1(), colind2_rand2());

size_type new_num_entries = rowptr2[rowptr2.size()-1];
TEST_EQUALITY(num_entries, new_num_entries);

TEST_COMPARE_ARRAYS(rowptr, rowptr2);

colind2_rand1.resize(num_entries);
TEST_COMPARE_ARRAYS(colind, colind2_rand1);

colind2_rand2.resize(num_entries);
TEST_COMPARE_ARRAYS(colind, colind2_rand2);

vals2.resize(num_entries);
vals2_rand.resize(num_entries);
for (size_type i=0; i<num_entries; i++) vals2[i] = 2.*vals[i];
TEST_COMPARE_FLOATING_ARRAYS(vals2, vals2_rand, 1.e-12);

}


TEUCHOS_UNIT_TEST_TEMPLATE_2_DECL( Import_Util,LowCommunicationMakeColMapAndReindex, LO, GO) {
// Test the colmap...
RCP<const Comm<int> > Comm = getDefaultComm();
Expand Down Expand Up @@ -2567,6 +2689,7 @@ TEUCHOS_UNIT_TEST_TEMPLATE_2_DECL( Import_Util,GetTwoTransferOwnershipVector, LO
TEUCHOS_UNIT_TEST_TEMPLATE_3_INSTANT( CrsMatrixImportExport, doImport, LO, GO, SC ) \
TEUCHOS_UNIT_TEST_TEMPLATE_3_INSTANT( FusedImportExport, doImport, LO, GO, SC ) \
TEUCHOS_UNIT_TEST_TEMPLATE_3_INSTANT( Import_Util, UnpackAndCombineWithOwningPIDs, LO, GO, SC ) \
TEUCHOS_UNIT_TEST_TEMPLATE_3_INSTANT( Import_Util, SortCrsEntries, LO, GO, SC ) \
TEUCHOS_UNIT_TEST_TEMPLATE_3_INSTANT( FusedImportExport, MueLuStyle, LO, GO, SC )

// Note: This test fails. Should fix later.
Expand Down

0 comments on commit d482822

Please sign in to comment.