Shuenhoy/vtk-m2
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'master' into particle_density

Browse Source
This commit is contained in:
Li-Ta Lo 2020-10-15 07:41:05 -06:00
commit 8c15745c04
29 changed files with 624 additions and 258 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

35
benchmarking/BenchmarkFilters.cxx
View File

@ -25,6 +25,7 @@
#include <vtkm/cont/Logging.h>
#include <vtkm/cont/RuntimeDeviceTracker.h>
#include <vtkm/cont/Timer.h>
#include <vtkm/cont/testing/MakeTestDataSet.h>
#include <vtkm/cont/internal/OptionParser.h>
@ -38,6 +39,7 @@
#include <vtkm/filter/Tetrahedralize.h>
#include <vtkm/filter/Threshold.h>
#include <vtkm/filter/ThresholdPoints.h>
#include <vtkm/filter/Triangulate.h>
#include <vtkm/filter/VectorMagnitude.h>
#include <vtkm/filter/VertexClustering.h>
#include <vtkm/filter/WarpScalar.h>
@ -91,12 +93,15 @@ vtkm::cont::InitializeResult Config;
// The input dataset we'll use on the filters:
static vtkm::cont::DataSet InputDataSet;
static vtkm::cont::DataSet UnstructuredInputDataSet;
// The point scalars to use:
static std::string PointScalarsName;
// The cell scalars to use:
static std::string CellScalarsName;
// The point vectors to use:
static std::string PointVectorsName;
// Whether the input is a file or is generated
bool FileAsInput = false;
bool InputIsStructured()
{
@ -398,7 +403,9 @@ void BenchContourGenerator(::benchmark::internal::Benchmark* bm)
helper(3);
helper(12);
}
VTKM_BENCHMARK_APPLY(BenchContour, BenchContourGenerator);
// :TODO: Disabled until SIGSEGV in Countour when passings field is resolved
//VTKM_BENCHMARK_APPLY(BenchContour, BenchContourGenerator);
void BenchExternalFaces(::benchmark::State& state)
{
@ -427,10 +434,9 @@ void BenchTetrahedralize(::benchmark::State& state)
const vtkm::cont::DeviceAdapterId device = Config.Device;
// This filter only supports structured datasets:
if (!InputIsStructured())
if (FileAsInput && !InputIsStructured())
{
state.SkipWithError("Tetrahedralize Filter requires structured data.");
return;
}
vtkm::filter::Tetrahedralize filter;
@ -455,10 +461,9 @@ void BenchVertexClustering(::benchmark::State& state)
const vtkm::Id numDivs = static_cast<vtkm::Id>(state.range(0));
// This filter only supports unstructured datasets:
if (InputIsStructured())
if (FileAsInput && InputIsStructured())
{
state.SkipWithError("VertexClustering Filter requires unstructured data.");
return;
state.SkipWithError("VertexClustering Filter requires unstructured data (use --tetra).");
}
vtkm::filter::VertexClustering filter;
@ -468,8 +473,9 @@ void BenchVertexClustering(::benchmark::State& state)
for (auto _ : state)
{
(void)_;
timer.Start();
auto result = filter.Execute(InputDataSet);
auto result = filter.Execute(UnstructuredInputDataSet);
::benchmark::DoNotOptimize(result);
timer.Stop();
@ -529,13 +535,12 @@ struct PrepareForInput
void BenchReverseConnectivityGen(::benchmark::State& state)
{
if (InputIsStructured())
if (FileAsInput && InputIsStructured())
{
state.SkipWithError("ReverseConnectivityGen requires unstructured data.");
return;
state.SkipWithError("ReverseConnectivityGen requires unstructured data (--use tetra).");
}
auto cellset = InputDataSet.GetCellSet();
auto cellset = UnstructuredInputDataSet.GetCellSet();
PrepareForInput functor;
for (auto _ : state)
{
@ -978,6 +983,7 @@ void InitDataSet(int& argc, char** argv)
std::cerr << "[InitDataSet] Loading file: " << filename << "\n";
vtkm::io::VTKDataSetReader reader(filename);
InputDataSet = reader.ReadDataSet();
FileAsInput = true;
}
else
{
@ -987,8 +993,15 @@ void InitDataSet(int& argc, char** argv)
source.SetExtent({ 0 }, { waveletDim - 1 });
InputDataSet = source.Execute();
vtkm::cont::DataSet input = vtkm::cont::testing::MakeTestDataSet().Make2DUniformDataSet2();
vtkm::filter::Triangulate triangulateFilter;
triangulateFilter.SetFieldsToPass(
vtkm::filter::FieldSelection(vtkm::filter::FieldSelection::MODE_ALL));
UnstructuredInputDataSet = triangulateFilter.Execute(input);
}
if (tetra)
{
std::cerr << "[InitDataSet] Tetrahedralizing dataset...\n";

3
data/data/third_party/ecl_cc/README vendored Normal file
View File

@ -0,0 +1,3 @@
version https://git-lfs.github.com/spec/v1
oid sha256:6f5e6e3dc559fefc7990daaec071fcd620f620e5ab8652dddaa6b43ca4ba08e7
size 222

3
data/data/third_party/ecl_cc/internet.egr vendored Normal file
View File

@ -0,0 +1,3 @@
version https://git-lfs.github.com/spec/v1
oid sha256:892470e152ccd46ddcca70e26bcd88816c247f08c496319cea80864b6b94ce46
size 3596536

4
vtkm/cont/ArrayHandleGroupVecVariable.h
View File

@ -459,6 +459,8 @@ VTKM_CONT void ConvertNumComponentsToOffsets(
using namespace vtkm::cont;
VTKM_IS_ARRAY_HANDLE(NumComponentsArrayType);
VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf);
Algorithm::ScanExtended(device, make_ArrayHandleCast<vtkm::Id>(numComponentsArray), offsetsArray);
sourceArraySize = ArrayGetValue(offsetsArray.GetNumberOfValues() - 1, offsetsArray);
@ -472,6 +474,8 @@ VTKM_CONT void ConvertNumComponentsToOffsets(
{
VTKM_IS_ARRAY_HANDLE(NumComponentsArrayType);
VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf);
vtkm::cont::Algorithm::ScanExtended(
device, vtkm::cont::make_ArrayHandleCast<vtkm::Id>(numComponentsArray), offsetsArray);
}

2
vtkm/cont/ArrayRangeCompute.hxx
View File

@ -47,6 +47,8 @@ inline vtkm::cont::ArrayHandle<vtkm::Range> ArrayRangeComputeImpl(
const vtkm::cont::ArrayHandle<T, S>& input,
vtkm::cont::DeviceAdapterId device)
{
VTKM_LOG_SCOPE(vtkm::cont::LogLevel::Perf, "ArrayRangeCompute");
using VecTraits = vtkm::VecTraits<T>;
using CT = typename VecTraits::ComponentType;
//We want to minimize the amount of code that we do in try execute as

2
vtkm/cont/CellLocatorBoundingIntervalHierarchy.cxx
View File

@ -229,6 +229,8 @@ CellLocatorBoundingIntervalHierarchy::~CellLocatorBoundingIntervalHierarchy() =
void CellLocatorBoundingIntervalHierarchy::Build()
{
VTKM_LOG_SCOPE(vtkm::cont::LogLevel::Perf, "CellLocatorBoundingIntervalHierarchy::Build");
vtkm::cont::Invoker invoker;
vtkm::cont::DynamicCellSet cellSet = this->GetCellSet();

2
vtkm/cont/CellLocatorTwoLevel.cxx
View File

@ -347,6 +347,8 @@ namespace cont
///
VTKM_CONT void CellLocatorTwoLevel::Build()
{
VTKM_LOG_SCOPE(vtkm::cont::LogLevel::Perf, "CellLocatorTwoLevel::Build");
vtkm::cont::Invoker invoke;
auto cellset = this->GetCellSet();

2
vtkm/cont/CellSetExplicit.h
View File

@ -81,6 +81,8 @@ template <typename S1, typename S2>
void ConvertNumIndicesToOffsets(const vtkm::cont::ArrayHandle<vtkm::Id, S1>& numIndices,
vtkm::cont::ArrayHandle<vtkm::Id, S2>& offsets)
{
VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf);
vtkm::cont::Algorithm::ScanExtended(numIndices, offsets);
}

2
vtkm/cont/Field.h
View File

@ -156,6 +156,8 @@ private:
{
VTKM_IS_LIST(TypeList);
VTKM_LOG_SCOPE(vtkm::cont::LogLevel::Perf, "Field::GetRange");
if (this->ModifiedFlag)
{
vtkm::cont::CastAndCall(

2
vtkm/cont/PointLocatorSparseGrid.cxx
View File

@ -61,6 +61,8 @@ private:
void PointLocatorSparseGrid::Build()
{
VTKM_LOG_SCOPE(vtkm::cont::LogLevel::Perf, "PointLocatorSparseGrid::Build");
if (this->IsRangeInvalid())
{
this->Range = this->GetCoordinates().GetRange();

11
vtkm/cont/internal/FunctorsGeneral.h
View File

@ -12,6 +12,7 @@
#include <vtkm/Algorithms.h>
#include <vtkm/BinaryOperators.h>
#include <vtkm/BinaryPredicates.h>
#include <vtkm/TypeTraits.h>
#include <vtkm/UnaryPredicates.h>
#include <vtkm/cont/ArrayPortalToIterators.h>
@ -83,15 +84,7 @@ struct WrappedBinaryOperator
}
};
//needs to be in a location that TBB DeviceAdapterAlgorithm can reach
struct DefaultCompareFunctor
{
template <typename T>
VTKM_EXEC bool operator()(const T& first, const T& second) const
{
return first < second;
}
};
using DefaultCompareFunctor = vtkm::SortLess;
//needs to be in a location that TBB DeviceAdapterAlgorithm can reach
template <typename T, typename U, class BinaryCompare = DefaultCompareFunctor>

54
vtkm/cont/kokkos/internal/DeviceAdapterAlgorithmKokkos.h
View File

@ -24,8 +24,11 @@
VTKM_THIRDPARTY_PRE_INCLUDE
#include <Kokkos_Core.hpp>
#include <Kokkos_DualView.hpp>
#include <Kokkos_Sort.hpp>
VTKM_THIRDPARTY_POST_INCLUDE
#include <type_traits>
namespace vtkm
{
namespace cont
@ -88,6 +91,10 @@ struct DeviceAdapterAlgorithm<vtkm::cont::DeviceAdapterTagKokkos>
vtkm::cont::DeviceAdapterTagKokkos>
{
private:
using Superclass = vtkm::cont::internal::DeviceAdapterAlgorithmGeneral<
DeviceAdapterAlgorithm<vtkm::cont::DeviceAdapterTagKokkos>,
vtkm::cont::DeviceAdapterTagKokkos>;
constexpr static vtkm::Id ErrorMessageMaxLength = 1024;
using ErrorMessageStorage =
Kokkos::DualView<char*, Kokkos::LayoutLeft, Kokkos::DefaultExecutionSpace>;
@ -134,6 +141,25 @@ public:
vtkm::Id{ 0 });
}
using Superclass::Copy;
template <typename T>
VTKM_CONT static void Copy(const vtkm::cont::ArrayHandle<T>& input,
vtkm::cont::ArrayHandle<T>& output)
{
const vtkm::Id inSize = input.GetNumberOfValues();
vtkm::cont::Token token;
auto portalIn = input.PrepareForInput(vtkm::cont::DeviceAdapterTagKokkos{}, token);
auto portalOut = output.PrepareForOutput(inSize, vtkm::cont::DeviceAdapterTagKokkos{}, token);
kokkos::internal::KokkosViewConstExec<T> viewIn(portalIn.GetArray(), inSize);
kokkos::internal::KokkosViewExec<T> viewOut(portalOut.GetArray(), inSize);
Kokkos::deep_copy(Kokkos::DefaultExecutionSpace{}, viewOut, viewIn);
}
template <typename WType, typename IType>
VTKM_CONT static void ScheduleTask(
vtkm::exec::kokkos::internal::TaskBasic1D<WType, IType>& functor,
@ -202,6 +228,34 @@ public:
ScheduleTask(kernel, rangeMax);
}
private:
template <typename T>
VTKM_CONT static void SortImpl(vtkm::cont::ArrayHandle<T>& values, vtkm::SortLess, std::true_type)
{
vtkm::cont::Token token;
auto portal = values.PrepareForInPlace(vtkm::cont::DeviceAdapterTagKokkos{}, token);
kokkos::internal::KokkosViewExec<T> view(portal.GetArray(), portal.GetNumberOfValues());
Kokkos::sort(view);
}
template <typename T>
VTKM_CONT static void SortImpl(vtkm::cont::ArrayHandle<T>& values,
vtkm::SortLess comp,
std::false_type)
{
Superclass::Sort(values, comp);
}
public:
using Superclass::Sort;
template <typename T>
VTKM_CONT static void Sort(vtkm::cont::ArrayHandle<T>& values, vtkm::SortLess comp)
{
SortImpl(values, comp, typename std::is_scalar<T>::type{});
}
VTKM_CONT static void Synchronize() {}
};

32
vtkm/cont/testing/TestingDeviceAdapter.h
View File

@ -81,25 +81,6 @@ private:
public:
// Cuda kernels have to be public (in Cuda 4.0).
struct CopyArrayKernel
{
VTKM_CONT
CopyArrayKernel(const IdPortalConstType& input, const IdPortalType& output)
: InputArray(input)
, OutputArray(output)
{
}
VTKM_EXEC void operator()(vtkm::Id index, const vtkm::exec::internal::ErrorMessageBuffer&) const
{
this->OutputArray.Set(index, this->InputArray.Get(index));
}
VTKM_CONT void SetErrorMessageBuffer(const vtkm::exec::internal::ErrorMessageBuffer&) {}
IdPortalConstType InputArray;
IdPortalType OutputArray;
};
template <typename PortalType>
struct GenericClearArrayKernel
@ -143,19 +124,6 @@ public:
using ClearArrayKernel = GenericClearArrayKernel<IdPortalType>;
struct ClearArrayMapKernel //: public vtkm::exec::WorkletMapField
{
// using ControlSignature = void(Field(Out));
// using ExecutionSignature = void(_1);
template <typename T>
VTKM_EXEC void operator()(T& value) const
{
value = OFFSET;
}
};
template <typename PortalType>
struct AddArrayKernel
{

4
vtkm/filter/MapFieldMergeAverage.cxx
View File

@ -40,6 +40,8 @@ bool vtkm::filter::MapFieldMergeAverage(const vtkm::cont::Field& inputField,
const vtkm::worklet::internal::KeysBase& keys,
vtkm::cont::Field& outputField)
{
VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf);
vtkm::cont::VariantArrayHandle outputArray;
DoMapFieldMerge functor;
inputField.GetData().ResetTypes<vtkm::TypeListAll>().CastAndCall(
@ -59,6 +61,8 @@ bool vtkm::filter::MapFieldMergeAverage(const vtkm::cont::Field& inputField,
const vtkm::worklet::internal::KeysBase& keys,
vtkm::cont::DataSet& outputData)
{
VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf);
vtkm::cont::Field outputField;
bool success = vtkm::filter::MapFieldMergeAverage(inputField, keys, outputField);
if (success)

2
vtkm/filter/MapFieldPermutation.cxx
View File

@ -78,6 +78,8 @@ VTKM_FILTER_COMMON_EXPORT VTKM_CONT bool vtkm::filter::MapFieldPermutation(
vtkm::cont::Field& outputField,
vtkm::Float64 invalidValue)
{
VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf);
vtkm::cont::VariantArrayHandle outputArray;
DoMapFieldPermutation functor;
inputField.GetData().ResetTypes<vtkm::TypeListAll>().CastAndCall(

13
vtkm/filter/Triangulate.hxx
View File

@ -17,13 +17,14 @@
namespace
{
class DeduceCellSet
class DeduceCellSetTriangulate
{
vtkm::worklet::Triangulate& Worklet;
vtkm::cont::CellSetSingleType<>& OutCellSet;
public:
DeduceCellSet(vtkm::worklet::Triangulate& worklet, vtkm::cont::CellSetSingleType<>& outCellSet)
DeduceCellSetTriangulate(vtkm::worklet::Triangulate& worklet,
vtkm::cont::CellSetSingleType<>& outCellSet)
: Worklet(worklet)
, OutCellSet(outCellSet)
{
@ -35,17 +36,17 @@ public:
}
};
template <>
void DeduceCellSet::operator()(const vtkm::cont::CellSetExplicit<>& cellset) const
void DeduceCellSetTriangulate::operator()(const vtkm::cont::CellSetExplicit<>& cellset) const
{
this->OutCellSet = Worklet.Run(cellset);
}
template <>
void DeduceCellSet::operator()(const vtkm::cont::CellSetStructured<2>& cellset) const
void DeduceCellSetTriangulate::operator()(const vtkm::cont::CellSetStructured<2>& cellset) const
{
this->OutCellSet = Worklet.Run(cellset);
}
template <>
void DeduceCellSet::operator()(const vtkm::cont::CellSetStructured<3>& cellset) const
void DeduceCellSetTriangulate::operator()(const vtkm::cont::CellSetStructured<3>& cellset) const
{
this->OutCellSet = Worklet.Run(cellset);
}
@ -72,7 +73,7 @@ inline VTKM_CONT vtkm::cont::DataSet Triangulate::DoExecute(
const vtkm::cont::DynamicCellSet& cells = input.GetCellSet();
vtkm::cont::CellSetSingleType<> outCellSet;
DeduceCellSet triangulate(this->Worklet, outCellSet);
DeduceCellSetTriangulate triangulate(this->Worklet, outCellSet);
vtkm::cont::CastAndCall(vtkm::filter::ApplyPolicyCellSet(cells, policy, *this), triangulate);

74
vtkm/filter/testing/UnitTestParticleAdvectionFilter.cxx
View File

@ -10,6 +10,7 @@
#include <vtkm/cont/DataSetBuilderUniform.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/filter/CleanGrid.h>
#include <vtkm/filter/ParticleAdvection.h>
#include <vtkm/io/VTKDataSetReader.h>
#include <vtkm/thirdparty/diy/environment.h>
@ -37,6 +38,8 @@ vtkm::cont::DataSet CreateDataSet(const vtkm::Id3& dims,
void TestBasic()
{
std::cout << "Basic uniform grid" << std::endl;
const vtkm::Id3 dims(5, 5, 5);
const vtkm::Vec3f origin(0, 0, 0), spacing(1, 1, 1), vecX(1, 0, 0);
vtkm::cont::DataSet ds = CreateDataSet(dims, origin, spacing, vecX);
@ -68,6 +71,8 @@ void TestBasic()
void TestPartitionedDataSet()
{
std::cout << "Partitioned data set" << std::endl;
const vtkm::Id3 dims(5, 5, 5);
const vtkm::Vec3f o1(0, 0, 0), o2(4, 0, 0), o3(8, 0, 0);
const vtkm::Vec3f spacing(1, 1, 1);
@ -95,13 +100,15 @@ void TestPartitionedDataSet()
particleAdvection.SetActiveField("vector");
auto out = particleAdvection.Execute(pds);
std::cout << "###### " << out.GetNumberOfPartitions() << std::endl;
#if 0
std::cout << "###### " << out.GetNumberOfPartitions() << std::endl;
for (int i = 0; i < out.GetNumberOfPartitions(); i++)
{
std::cout << "PID= " << i << std::endl;
out.GetPartition(i).PrintSummary(std::cout);
}
#endif
VTKM_TEST_ASSERT(out.GetNumberOfPartitions() == 1, "Wrong number of partitions in output");
auto ds = out.GetPartition(0);
@ -119,30 +126,15 @@ void TestPartitionedDataSet()
vtkm::cont::DynamicCellSet dcells = ds.GetCellSet();
VTKM_TEST_ASSERT(dcells.GetNumberOfCells() == numSeeds, "Wrong number of cells");
ds.PrintSummary(std::cout);
//ds.PrintSummary(std::cout);
}
void TestFile(const std::string& fname,
const std::vector<vtkm::Vec3f>& pts,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps,
const std::vector<vtkm::Vec3f>& endPts)
void TestDataSet(const vtkm::cont::DataSet& dataset,
const std::vector<vtkm::Vec3f>& pts,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps,
const std::vector<vtkm::Vec3f>& endPts)
{
vtkm::io::VTKDataSetReader reader(fname);
vtkm::cont::DataSet ds;
try
{
ds = reader.ReadDataSet();
}
catch (vtkm::io::ErrorIO& e)
{
std::string message("Error reading: ");
message += fname;
message += ", ";
message += e.GetMessage();
VTKM_TEST_FAIL(message.c_str());
}
vtkm::Id numPoints = static_cast<vtkm::Id>(pts.size());
std::vector<vtkm::Particle> seeds;
@ -156,7 +148,7 @@ void TestFile(const std::string& fname,
particleAdvection.SetSeeds(seedArray);
particleAdvection.SetActiveField("vec");
auto output = particleAdvection.Execute(ds);
auto output = particleAdvection.Execute(dataset);
auto coords = output.GetCoordinateSystem().GetDataAsMultiplexer();
vtkm::cont::DynamicCellSet dcells = output.GetCellSet();
@ -174,6 +166,42 @@ void TestFile(const std::string& fname,
}
}
void TestFile(const std::string& fname,
const std::vector<vtkm::Vec3f>& pts,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps,
const std::vector<vtkm::Vec3f>& endPts)
{
std::cout << fname << std::endl;
vtkm::io::VTKDataSetReader reader(fname);
vtkm::cont::DataSet dataset;
try
{
dataset = reader.ReadDataSet();
}
catch (vtkm::io::ErrorIO& e)
{
std::string message("Error reading: ");
message += fname;
message += ", ";
message += e.GetMessage();
VTKM_TEST_FAIL(message.c_str());
}
TestDataSet(dataset, pts, stepSize, maxSteps, endPts);
std::cout << " as explicit grid" << std::endl;
vtkm::filter::CleanGrid clean;
clean.SetCompactPointFields(false);
clean.SetMergePoints(false);
clean.SetRemoveDegenerateCells(false);
vtkm::cont::DataSet explicitData = clean.Execute(dataset);
TestDataSet(explicitData, pts, stepSize, maxSteps, endPts);
}
void TestParticleAdvectionFilter()
{
TestBasic();

3
vtkm/worklet/AverageByKey.h
View File

@ -70,6 +70,7 @@ struct AverageByKey
const InArrayType& inValues,
OutArrayType& outAverages)
{
VTKM_LOG_SCOPE(vtkm::cont::LogLevel::Perf, "AverageByKey::Run");
vtkm::worklet::DispatcherReduceByKey<AverageWorklet> dispatcher;
dispatcher.Invoke(keys, inValues, outAverages);
@ -133,6 +134,8 @@ struct AverageByKey
vtkm::cont::ArrayHandle<KeyType, KeyOutStorage>& outputKeyArray,
vtkm::cont::ArrayHandle<ValueType, ValueOutStorage>& outputValueArray)
{
VTKM_LOG_SCOPE(vtkm::cont::LogLevel::Perf, "AverageByKey::Run");
auto results = vtkm::worklet::DescriptiveStatistics::Run(keyArray, valueArray);
// Copy/TransformCopy from results to outputKeyArray and outputValueArray

2
vtkm/worklet/MaskSelect.cxx
View File

@ -130,6 +130,8 @@ vtkm::worklet::MaskSelect::ThreadToOutputMapType vtkm::worklet::MaskSelect::Buil
const VariantArrayHandleMask& maskArray,
vtkm::cont::DeviceAdapterId device)
{
VTKM_LOG_SCOPE(vtkm::cont::LogLevel::Perf, "MaskSelect::Build");
vtkm::worklet::MaskSelect::ThreadToOutputMapType threadToOutputMap;
maskArray.CastAndCall(MaskBuilder(), threadToOutputMap, device);
return threadToOutputMap;

2
vtkm/worklet/ScatterCounting.cxx
View File

@ -203,6 +203,8 @@ void vtkm::worklet::ScatterCounting::BuildArrays(const VariantArrayHandleCount&
vtkm::cont::DeviceAdapterId device,
bool saveInputToOutputMap)
{
VTKM_LOG_SCOPE(vtkm::cont::LogLevel::Perf, "ScatterCounting::BuildArrays");
countArray.CastAndCall(
vtkm::worklet::detail::ScatterCountingBuilder(), device, saveInputToOutputMap, this);
}

96
vtkm/worklet/connectivities/GraphConnectivity.h
View File

@ -24,94 +24,74 @@ namespace connectivity
{
namespace detail
{
class Graft : public vtkm::worklet::WorkletMapField
class GraphGraft : public vtkm::worklet::WorkletMapField
{
public:
using ControlSignature = void(FieldIn start,
FieldIn degree,
WholeArrayIn ids,
WholeArrayInOut comp);
AtomicArrayInOut comp);
using ExecutionSignature = void(WorkIndex, _1, _2, _3, _4);
using InputDomain = _1;
// TODO: Use Scatter?
template <typename InPortalType, typename InOutPortalType>
template <typename InPortalType, typename AtomicCompInOut>
VTKM_EXEC void operator()(vtkm::Id index,
vtkm::Id start,
vtkm::Id degree,
const InPortalType& conn,
InOutPortalType& comp) const
AtomicCompInOut& comp) const
{
for (vtkm::Id offset = start; offset < start + degree; offset++)
{
vtkm::Id neighbor = conn.Get(offset);
if ((comp.Get(index) == comp.Get(comp.Get(index))) && (comp.Get(neighbor) < comp.Get(index)))
{
comp.Set(comp.Get(index), comp.Get(neighbor));
}
// We need to reload thisComp and thatComp every iteration since
// they might have been changed by Unite() both as a result of
// attaching one tree to the other or as a result of path compression
// in findRoot().
auto thisComp = comp.Get(index);
auto thatComp = comp.Get(neighbor);
// Merge the two components one way or the other, the order will
// be resolved by Unite().
UnionFind::Unite(comp, thisComp, thatComp);
}
}
};
}
// Single pass connected component algorithm from
// Jaiganesh, Jayadharini, and Martin Burtscher.
// "A high-performance connected components implementation for GPUs."
// Proceedings of the 27th International Symposium on High-Performance
// Parallel and Distributed Computing. 2018.
class GraphConnectivity
{
public:
using Algorithm = vtkm::cont::Algorithm;
template <typename InputPortalType, typename OutputPortalType>
void Run(const InputPortalType& numIndicesArray,
const InputPortalType& indexOffsetsArray,
const InputPortalType& connectivityArray,
OutputPortalType& componentsOut) const
template <typename InputArrayType, typename OutputArrayType>
void Run(const InputArrayType& numIndicesArray,
const InputArrayType& indexOffsetsArray,
const InputArrayType& connectivityArray,
OutputArrayType& componentsOut) const
{
bool everythingIsAStar = false;
vtkm::cont::ArrayHandle<vtkm::Id> components;
Algorithm::Copy(
vtkm::cont::ArrayHandleCounting<vtkm::Id>(0, 1, numIndicesArray.GetNumberOfValues()),
components);
VTKM_IS_ARRAY_HANDLE(InputArrayType);
VTKM_IS_ARRAY_HANDLE(OutputArrayType);
//used as an atomic bool, so we use Int32 as it the
//smallest type that VTK-m supports as atomics
vtkm::cont::ArrayHandle<vtkm::Int32> allStars;
allStars.Allocate(1);
using Algorithm = vtkm::cont::Algorithm;
// Initialize the parent pointer to point to the node itself. There are other
// ways to initialize the parent pointers, for example, a smaller or the minimal
// neighbor.
Algorithm::Copy(vtkm::cont::ArrayHandleIndex(numIndicesArray.GetNumberOfValues()),
componentsOut);
vtkm::cont::Invoker invoke;
do
{
allStars.WritePortal().Set(0, 1); //reset the atomic state
invoke(detail::Graft{}, indexOffsetsArray, numIndicesArray, connectivityArray, components);
// Detection of allStars has to come before pointer jumping. Don't try to rearrange it.
invoke(IsStar{}, components, allStars);
everythingIsAStar = (allStars.WritePortal().Get(0) == 1);
invoke(PointerJumping{}, components);
} while (!everythingIsAStar);
invoke(
detail::GraphGraft{}, indexOffsetsArray, numIndicesArray, connectivityArray, componentsOut);
invoke(PointerJumping{}, componentsOut);
// renumber connected component to the range of [0, number of components).
vtkm::cont::ArrayHandle<vtkm::Id> uniqueComponents;
Algorithm::Copy(components, uniqueComponents);
Algorithm::Sort(uniqueComponents);
Algorithm::Unique(uniqueComponents);
vtkm::cont::ArrayHandle<vtkm::Id> cellIds;
Algorithm::Copy(
vtkm::cont::ArrayHandleCounting<vtkm::Id>(0, 1, numIndicesArray.GetNumberOfValues()),
cellIds);
vtkm::cont::ArrayHandle<vtkm::Id> uniqueColor;
Algorithm::Copy(
vtkm::cont::ArrayHandleCounting<vtkm::Id>(0, 1, uniqueComponents.GetNumberOfValues()),
uniqueColor);
vtkm::cont::ArrayHandle<vtkm::Id> cellColors;
vtkm::cont::ArrayHandle<vtkm::Id> cellIdsOut;
InnerJoin().Run(
components, cellIds, uniqueComponents, uniqueColor, cellColors, cellIdsOut, componentsOut);
Algorithm::SortByKey(cellIdsOut, componentsOut);
Renumber::Run(componentsOut);
}
};
}

118
vtkm/worklet/connectivities/ImageConnectivity.h
View File

@ -32,69 +32,59 @@ class ImageGraft : public vtkm::worklet::WorkletPointNeighborhood
{
public:
using ControlSignature = void(CellSetIn,
FieldInNeighborhood compIn,
FieldInNeighborhood color,
WholeArrayInOut compOut,
AtomicArrayInOut changed);
FieldInNeighborhood neighborComp,
FieldInNeighborhood neighborColor,
AtomicArrayInOut compOut);
using ExecutionSignature = void(WorkIndex, _2, _3, _4, _5);
using ExecutionSignature = void(Boundary, _2, _3, _4);
template <typename Comp>
VTKM_EXEC vtkm::Id findRoot(Comp& comp, vtkm::Id index) const
{
while (comp.Get(index) != index)
index = comp.Get(index);
return index;
}
// compOut is an alias of compIn such that we can update component labels
template <typename NeighborComp, typename NeighborColor, typename CompOut, typename AtomicInOut>
VTKM_EXEC void operator()(const vtkm::Id index,
// compOut is a "linear" alias of neightborComp such that we can update component labels
template <typename Boundary,
typename NeighborComp,
typename NeighborColor,
typename AtomicCompOut>
VTKM_EXEC void operator()(Boundary boundary,
const NeighborComp& neighborComp,
const NeighborColor& neighborColor,
CompOut& compOut,
AtomicInOut& updated) const
AtomicCompOut& compOut) const
{
vtkm::Id myComp = neighborComp.Get(0, 0, 0);
auto minComp = myComp;
auto thisColor = neighborColor.Get(0, 0, 0);
auto myColor = neighborColor.Get(0, 0, 0);
auto minIndices = boundary.MinNeighborIndices(1);
auto maxIndices = boundary.MaxNeighborIndices(1);
for (int k = -1; k <= 1; k++)
for (int k = minIndices[2]; k <= maxIndices[2]; k++)
{
for (int j = -1; j <= 1; j++)
for (int j = minIndices[1]; j <= maxIndices[1]; j++)
{
for (int i = -1; i <= 1; i++)
for (int i = minIndices[0]; i <= maxIndices[0]; i++)
{
if (myColor == neighborColor.Get(i, j, k))
if (thisColor == neighborColor.Get(i, j, k))
{
minComp = vtkm::Min(minComp, neighborComp.Get(i, j, k));
// We need to reload thisComp and thatComp every iteration since
// they might have been changed by Unite(), both as a result of
// attaching one tree to the other or as a result of path compaction
// in findRoot().
auto thisComp = neighborComp.Get(0, 0, 0);
auto thatComp = neighborComp.Get(i, j, k);
// Merge the two components one way or the other, the order will
// be resolved by Unite().
UnionFind::Unite(compOut, thisComp, thatComp);
}
}
}
}
// I don't just only want to update the component label of this pixel, I actually
// want to Union(FindRoot(myComponent), FindRoot(minComp)) and then Flatten the
// result.
compOut.Set(index, minComp);
auto myRoot = findRoot(compOut, myComp);
auto newRoot = findRoot(compOut, minComp);
if (myRoot < newRoot)
compOut.Set(newRoot, myRoot);
else if (myRoot > newRoot)
compOut.Set(myRoot, newRoot);
// mark an update occurred if no other worklets have done so yet
if (myComp != minComp && updated.Get(0) == 0)
{
updated.Set(0, 1);
}
}
};
}
// Single pass connected component algorithm from
// Jaiganesh, Jayadharini, and Martin Burtscher.
// "A high-performance connected components implementation for GPUs."
// Proceedings of the 27th International Symposium on High-Performance
// Parallel and Distributed Computing. 2018.
class ImageConnectivity
{
public:
@ -104,47 +94,21 @@ public:
template <int Dimension, typename T, typename StorageT, typename OutputPortalType>
void operator()(const vtkm::cont::ArrayHandle<T, StorageT>& pixels,
const vtkm::cont::CellSetStructured<Dimension>& input,
OutputPortalType& components) const
OutputPortalType& componentsOut) const
{
using Algorithm = vtkm::cont::Algorithm;
Algorithm::Copy(vtkm::cont::ArrayHandleCounting<vtkm::Id>(0, 1, pixels.GetNumberOfValues()),
components);
//used as an atomic bool, so we use Int32 as it the
//smallest type that VTK-m supports as atomics
vtkm::cont::ArrayHandle<vtkm::Int32> updateRequired;
updateRequired.Allocate(1);
// Initialize the parent pointer to point to the pixel itself. There are other
// ways to initialize the parent pointers, for example, a smaller or the minimal
// neighbor.
Algorithm::Copy(vtkm::cont::ArrayHandleIndex(pixels.GetNumberOfValues()), componentsOut);
vtkm::cont::Invoker invoke;
do
{
updateRequired.WritePortal().Set(0, 0); //reset the atomic state
invoke(detail::ImageGraft{}, input, components, pixels, components, updateRequired);
invoke(PointerJumping{}, components);
} while (updateRequired.WritePortal().Get(0) > 0);
invoke(detail::ImageGraft{}, input, componentsOut, pixels, componentsOut);
invoke(PointerJumping{}, componentsOut);
// renumber connected component to the range of [0, number of components).
vtkm::cont::ArrayHandle<vtkm::Id> uniqueComponents;
Algorithm::Copy(components, uniqueComponents);
Algorithm::Sort(uniqueComponents);
Algorithm::Unique(uniqueComponents);
vtkm::cont::ArrayHandle<vtkm::Id> pixelIds;
Algorithm::Copy(vtkm::cont::ArrayHandleCounting<vtkm::Id>(0, 1, pixels.GetNumberOfValues()),
pixelIds);
vtkm::cont::ArrayHandle<vtkm::Id> uniqueColor;
Algorithm::Copy(
vtkm::cont::ArrayHandleCounting<vtkm::Id>(0, 1, uniqueComponents.GetNumberOfValues()),
uniqueColor);
vtkm::cont::ArrayHandle<vtkm::Id> cellColors;
vtkm::cont::ArrayHandle<vtkm::Id> pixelIdsOut;
InnerJoin().Run(
components, pixelIds, uniqueComponents, uniqueColor, cellColors, pixelIdsOut, components);
Algorithm::SortByKey(pixelIdsOut, components);
Renumber::Run(componentsOut);
}
};

53
vtkm/worklet/connectivities/InnerJoin.h
View File

@ -11,7 +11,7 @@
#ifndef vtk_m_worklet_connectivity_InnerJoin_h
#define vtk_m_worklet_connectivity_InnerJoin_h
#include <vtkm/cont/ArrayHandleCounting.h>
#include <vtkm/cont/ArrayHandleIndex.h>
#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/worklet/ScatterCounting.h>
#include <vtkm/worklet/WorkletMapField.h>
@ -56,13 +56,13 @@ public:
// TODO: not mutating input keys and values?
template <typename Key, typename Value1, typename Value2>
void Run(vtkm::cont::ArrayHandle<Key>& key1,
vtkm::cont::ArrayHandle<Value1>& value1,
vtkm::cont::ArrayHandle<Key>& key2,
vtkm::cont::ArrayHandle<Value2>& value2,
vtkm::cont::ArrayHandle<Key>& keyOut,
vtkm::cont::ArrayHandle<Value1>& value1Out,
vtkm::cont::ArrayHandle<Value2>& value2Out) const
static void Run(vtkm::cont::ArrayHandle<Key>& key1,
vtkm::cont::ArrayHandle<Value1>& value1,
vtkm::cont::ArrayHandle<Key>& key2,
vtkm::cont::ArrayHandle<Value2>& value2,
vtkm::cont::ArrayHandle<Key>& keyOut,
vtkm::cont::ArrayHandle<Value1>& value1Out,
vtkm::cont::ArrayHandle<Value2>& value2Out)
{
Algorithm::SortByKey(key1, value1);
Algorithm::SortByKey(key2, value2);
@ -80,6 +80,43 @@ public:
mergeDisp.Invoke(key1, value1, lbs, value2, keyOut, value1Out, value2Out);
}
};
class Renumber
{
public:
static void Run(vtkm::cont::ArrayHandle<vtkm::Id>& componentsInOut)
{
using Algorithm = vtkm::cont::Algorithm;
// FIXME: we should be able to apply findRoot to each pixel and use some kind
// of atomic operation to get the number of unique components without the
// cost of copying and sorting. This might be able to be extended to also
// work for the renumbering (replacing InnerJoin) through atomic increment.
vtkm::cont::ArrayHandle<vtkm::Id> uniqueComponents;
Algorithm::Copy(componentsInOut, uniqueComponents);
Algorithm::Sort(uniqueComponents);
Algorithm::Unique(uniqueComponents);
vtkm::cont::ArrayHandle<vtkm::Id> ids;
Algorithm::Copy(vtkm::cont::ArrayHandleIndex(componentsInOut.GetNumberOfValues()), ids);
vtkm::cont::ArrayHandle<vtkm::Id> uniqueColor;
Algorithm::Copy(vtkm::cont::ArrayHandleIndex(uniqueComponents.GetNumberOfValues()),
uniqueColor);
vtkm::cont::ArrayHandle<vtkm::Id> cellColors;
vtkm::cont::ArrayHandle<vtkm::Id> pixelIdsOut;
InnerJoin::Run(componentsInOut,
ids,
uniqueComponents,
uniqueColor,
cellColors,
pixelIdsOut,
componentsInOut);
Algorithm::SortByKey(pixelIdsOut, componentsInOut);
}
};
}
}
} // vtkm::worklet::connectivity

171
vtkm/worklet/connectivities/UnionFind.h
View File

@ -20,6 +20,143 @@ namespace worklet
namespace connectivity
{
// Reference:
// Jayanti, Siddhartha V., and Robert E. Tarjan.
// "Concurrent Disjoint Set Union." arXiv preprint arXiv:2003.01203 (2020).
class UnionFind
{
public:
// This is the naive findRoot() without path compaction in SV Jayanti et. al.
// Since the parents array is read-only in this function, there is no data
// race when it is called by multiple treads concurrently. We can just call
// Get(). For cases where findRoot() is used in other functions that write
// to parents (e.g. Unite()), Get() actually calls Load() with
// memory_order_acquire ordering, which in turn ensure writes by other threads
// are reflected.
template <typename Parents>
static VTKM_EXEC vtkm::Id findRoot(const Parents& parents, vtkm::Id index)
{
while (parents.Get(index) != index)
index = parents.Get(index);
return index;
}
template <typename Parents>
static VTKM_EXEC void Unite(Parents& parents, vtkm::Id u, vtkm::Id v)
{
// Data Race Resolutions
// Since this function modifies the Union-Find data structure, concurrent
// invocation of it by 2 or more threads causes potential data race. Here
// is a case analysis why the potential data race does no harm in the
// context of the single pass connected component algorithm.
// Case 1, Two threads calling Unite(u, v) (and/or Unite(v, u)) concurrently.
// Problem: One thread might attach u to v while the other thread attach
// v to u, causing a cycle in the Union-Find data structure.
// Resolution: This is not so much of a race condition but a problem with
// the consistency of the algorithm. This can also happen in serial.
// This is resolved by "linking by index" as in SV Jayanti et.al. with less
// than as the total order. The two threads will make the same decision on
// how to Unite the two trees (e.g. from root with larger id to root with
// smaller id.) This avoids cycles in the resulting graph and maintains the
// rooted forest structure of Union-Find at the expense of duplicated (but
// benign) work.
// Case 2, T0 calling Unite(u, v) and T1 calling Unite(u, w) concurrently.
// Problem I: There is a potential write after read data race. After T0
// calls findRoot() for u but before actually updating the parent of root_u,
// T1 might have changed root_u to root_w and made root_u "obsolete".
// When the root of the tree to be attached to (e.g. root_u, when root_u < root_v)
// is changed, there is no hazard, since we are just attaching a tree to a
// now a non-root node, root_u, (thus, root_w <- root_u <- root_v).
// However, when the root of the attaching tree (root_v) is changed, it
// means that the root_u has been attached to yet some other root_s and became
// a non-root node. If we are now attaching this non-root node to root_w we
// would leave root_s behind and undoing previous work.
// Atomic Load/Store with memory_order_acquire are not able to detect this
// data race. While Load sees all previous Stores by other threads, it can not
// be aware of any Store after the Load.
// Resolution: Use atomic Compare and Swap in a loop when updating root_u.
// CAS will check if root of u has been updated by some other thread between
// findRoot(u) is called and when root of u is going to be updated. This is
// done by comparing the root_u = findRoot(u) to the current value at
// parents[root_u]. If they are the same, no data race has happened and the
// value in parents[root_u] is updated. However, if root_u != parent[root_u],
// it means parent[root_u] has been updated by some other thread. CAS returns
// the new value of parent[root_u] (root_s in the Problem description) which
// we can use as the new root_u. We keep retrying until there is no more
// data race and root_u == root_v i.e. they are in the same component.
// Problem II: There is a potential concurrent write data race as it is
// possible for the two threads to try to change the same old root to
// different new roots, e.g. T0 calls parents.Set(root_u, root_v) while T1
// calls parents.Set(root_u, root_w) where root_v < root_u and root_w < root_u
// (but the order of root_v and root_w is unspecified.) Each thread assumes
// success while the outcome is actually unspecified.
// Resolution: Use an atomic Compare and Swap is suggested in SV Janati et. al.
// as well as J. Jaiganesht et. al. to resolve the data race. The CAS
// checks if the old root is the same as what we expected. If so, there is
// no data race, CAS will set the root to the desired new value. The return
// value from CAS will equal to our expected old root and signifies a
// successful write which terminates the while loop.
// If the old root is not what we expected, it has been updated by some
// other thread and the update by this thread fails. The root as updated by
// the other thread is returned. This returned value would not equal to
// our desired new root, signifying the need to retry with the while loop.
// We can use this return "new root" as is without calling findRoot() to
// find the "new root". The while loop terminates when both u and v have
// the same root (thus united).
auto root_u = UnionFind::findRoot(parents, u);
auto root_v = UnionFind::findRoot(parents, v);
while (root_u != root_v)
{
// FIXME: we might be executing the loop one extra time than necessary.
// Nota Bene: VTKm's CompareAndSwap has a different order of parameters
// than common practice, it is (index, new, expected) rather than
// (index, expected, new).
if (root_u < root_v)
root_v = parents.CompareAndSwap(root_v, root_u, root_v);
else if (root_u > root_v)
root_u = parents.CompareAndSwap(root_u, root_v, root_u);
}
}
// This compresses the path from each node to its root thus flattening the
// trees and guarantees that the output trees will be rooted stars, i.e.
// they all have depth of 1.
// There is a "seemly" data race for this function. The root returned by
// findRoot() in one thread might become out of date if some other
// thread changed it before this function calls parents.Set() making the
// result tree not "short" enough and thus calls for a CompareAndSwap retry
// loop. However, this data race does not happen for the following reasons:
// 1. Since the only way for a root of a tree to be changed is through Unite(),
// as long as there is no concurrent invocation of Unite() and Flatten() there
// is no data race. This applies even for a compacting findRoot() which can
// still only change the parents of non-root nodes.
// 2. By the same token, since findRoot() does not change root and most
// "damage" parents.Set() can do is resetting root's parent to itself,
// the root of a tree can never be changed by this function. Thus, here
// is no data race between concurrent invocations of this function.
// Since the current findRoot() does not do path compaction, this algorithm
// has O(n) depth with O(n^2) of total work on a Parallel Random Access
// Machine (PRAM). However, we don't live in a synchronous, infinite number
// of processor PRAM world. In reality, since we put "parent pointers" in a
// array and all the pointers are pointing from larger indices to smaller
// ones, invocation for nodes with smaller ids are mostly likely be scheduled
// before and completes earlier than nodes with larger ids. This makes
// the "effective" path length shorter for nodes with larger ids.
// In this way, concurrency actually helps with algorithm complexity.
template <typename Parents>
static VTKM_EXEC void Flatten(Parents& parents, vtkm::Id index)
{
auto root = findRoot(parents, index);
parents.Set(index, root);
}
};
class PointerJumping : public vtkm::worklet::WorkletMapField
{
public:
@ -27,40 +164,10 @@ public:
using ExecutionSignature = void(WorkIndex, _1);
using InputDomain = _1;
template <typename Comp>
VTKM_EXEC vtkm::Id findRoot(Comp& comp, vtkm::Id index) const
{
while (comp.Get(index) != index)
index = comp.Get(index);
return index;
}
template <typename InOutPortalType>
VTKM_EXEC void operator()(vtkm::Id index, InOutPortalType& comp) const
VTKM_EXEC void operator()(vtkm::Id index, InOutPortalType& comps) const
{
// TODO: is there a data race between findRoot and comp.Set?
auto root = findRoot(comp, index);
comp.Set(index, root);
}
};
class IsStar : public vtkm::worklet::WorkletMapField
{
public:
using ControlSignature = void(WholeArrayIn comp, AtomicArrayInOut);
using ExecutionSignature = void(WorkIndex, _1, _2);
using InputDomain = _1;
template <typename InOutPortalType, typename AtomicInOut>
VTKM_EXEC void operator()(vtkm::Id index, InOutPortalType& comp, AtomicInOut& hasStar) const
{
//hasStar emulates a LogicalAnd across all the values
//where we start with a value of 'true'|1.
const bool isAStar = (comp.Get(index) == comp.Get(comp.Get(index)));
if (!isAStar && hasStar.Get(0) == 1)
{
hasStar.Set(0, 0);
}
UnionFind::Flatten(comps, index);
}
};

2
vtkm/worklet/internal/DispatcherBase.h
View File

@ -599,7 +599,7 @@ public:
{
VTKM_LOG_SCOPE(vtkm::cont::LogLevel::Perf,
"Invoking Worklet: '%s'",
vtkm::cont::TypeToString<WorkletType>().c_str());
vtkm::cont::TypeToString<DerivedClass>().c_str());
this->StartInvoke(std::forward<Args>(args)...);
}

1
vtkm/worklet/testing/CMakeLists.txt
View File

@ -49,6 +49,7 @@ set(unit_tests
UnitTestOrientNormals.cxx
UnitTestParticleAdvection.cxx
UnitTestPointElevation.cxx
UnitTestPointerJumping.cxx
UnitTestPointGradient.cxx
UnitTestPointTransform.cxx
UnitTestProbe.cxx

103
vtkm/worklet/testing/UnitTestGraphConnectivity.cxx
View File

@ -12,10 +12,105 @@
#include <vtkm/worklet/connectivities/GraphConnectivity.h>
class AdjacentDifference : public vtkm::worklet::WorkletMapField
{
public:
using ControlSignature = void(FieldIn index, WholeArrayIn counts, FieldOut outputCount);
using ExecutionSignature = void(_1, _2, _3);
using InputDomain = _1;
template <typename WholeArrayType>
VTKM_EXEC void operator()(const vtkm::Id& index,
const WholeArrayType& counts,
int& difference) const
{
difference = counts.Get(index + 1) - counts.Get(index);
}
};
class SameComponent : public vtkm::worklet::WorkletMapField
{
public:
using ControlSignature = void(FieldIn start,
FieldIn degree,
WholeArrayIn conns,
WholeArrayIn comps,
AtomicArrayInOut same);
using ExecutionSignature = void(WorkIndex, _1, _2, _3, _4, _5);
template <typename Conn, typename Comp, typename AtomicSame>
VTKM_EXEC void operator()(vtkm::Id index,
int start,
int degree,
const Conn& conns,
const Comp& comps,
AtomicSame& same) const
{
for (vtkm::Id offset = start; offset < start + degree; ++offset)
{
vtkm::Id neighbor = conns.Get(offset);
if (comps.Get(index) != comps.Get(neighbor))
{
same.Set(0, 0);
}
}
}
};
class TestGraphConnectivity
{
public:
void operator()() const
void TestECL_CC(const std::string& filename, int ncomps) const
{
auto pathname =
vtkm::cont::testing::Testing::GetTestDataBasePath() + "/third_party/ecl_cc/" + filename;
std::ifstream stream(pathname, std::ios_base::in | std::ios_base::binary);
int nnodes;
stream.read(reinterpret_cast<char*>(&nnodes), sizeof(nnodes));
int nedges;
stream.read(reinterpret_cast<char*>(&nedges), sizeof(nedges));
// CSR, there is one more element in offsets thant the actual number of nodes.
std::vector<int> offsets(nnodes + 1);
std::vector<int> conns(nedges);
stream.read(reinterpret_cast<char*>(offsets.data()), (nnodes + 1) * sizeof(int));
stream.read(reinterpret_cast<char*>(conns.data()), nedges * sizeof(int));
vtkm::cont::ArrayHandle<int> counts_h;
vtkm::cont::Invoker invoke;
invoke(AdjacentDifference{},
vtkm::cont::make_ArrayHandleCounting(0, 1, nnodes),
vtkm::cont::make_ArrayHandle<int>(offsets, vtkm::CopyFlag::On),
counts_h);
offsets.pop_back();
vtkm::cont::ArrayHandle<int> offsets_h =
vtkm::cont::make_ArrayHandle(offsets, vtkm::CopyFlag::On);
vtkm::cont::ArrayHandle<int> conns_h = vtkm::cont::make_ArrayHandle(conns, vtkm::CopyFlag::Off);
vtkm::cont::ArrayHandle<vtkm::Id> comps_h;
vtkm::worklet::connectivity::GraphConnectivity().Run(counts_h, offsets_h, conns_h, comps_h);
VTKM_TEST_ASSERT(vtkm::cont::Algorithm::Reduce(comps_h, vtkm::Id(0), vtkm::Maximum{}) ==
ncomps - 1,
"number of components mismatch");
vtkm::cont::ArrayHandle<vtkm::UInt32> atomicSame;
atomicSame.Allocate(1);
atomicSame.WritePortal().Set(0, 1);
invoke(SameComponent{}, offsets_h, counts_h, conns_h, comps_h, atomicSame);
VTKM_TEST_ASSERT(atomicSame.ReadPortal().Get(0) == 1,
"Neighboring nodes don't have the same component id");
}
void TestECL_CC_DataSets() const { TestECL_CC("internet.egr", 1); }
void TestSimpleGraph() const
{
vtkm::cont::ArrayHandle<vtkm::Id> counts_h =
vtkm::cont::make_ArrayHandle<vtkm::Id>({ 1, 1, 2, 2, 2 });
@ -32,6 +127,12 @@ public:
VTKM_TEST_ASSERT(comps.ReadPortal().Get(i) == 0, "Components has unexpected value.");
}
}
void operator()() const
{
TestSimpleGraph();
TestECL_CC_DataSets();
}
};
int UnitTestGraphConnectivity(int argc, char* argv[])

45
vtkm/worklet/testing/UnitTestImageConnectivity.cxx
View File

@ -11,6 +11,7 @@
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/filter/Contour.h>
#include <iomanip>
#include <vtkm/worklet/connectivities/ImageConnectivity.h>
class TestImageConnectivity
@ -22,6 +23,7 @@ public:
{
CCL_CUDA8x4();
CCL_CUDA8x8();
Valentine();
}
void CCL_CUDA8x4() const
@ -88,6 +90,49 @@ public:
"Components has unexpected value.");
}
}
void Valentine() const
{
// Sample image by VALENTINE PELTIER
// clang-format off
auto pixels = vtkm::cont::make_ArrayHandle<vtkm::UInt8>( {
1, 1, 0, 1, 0, 0,
0, 0, 0, 1, 1, 0,
1, 1, 0, 1, 0, 1,
1, 0, 1, 0, 0, 0,
0, 1, 0, 1, 1, 1,
1, 1, 0, 0, 1, 0,
});
// clang-format on
vtkm::cont::DataSetBuilderUniform builder;
vtkm::cont::DataSet data = builder.Create(vtkm::Id3(6, 6, 1));
auto colorField = vtkm::cont::make_FieldPoint("color", pixels);
data.AddField(colorField);
vtkm::cont::ArrayHandle<vtkm::Id> component;
vtkm::worklet::connectivity::ImageConnectivity().Run(
data.GetCellSet().Cast<vtkm::cont::CellSetStructured<2>>(), colorField.GetData(), component);
// clang-format off
std::vector<vtkm::UInt8> componentExpected = {
0, 0, 1, 2, 1, 1,
1, 1, 1, 2, 2, 1,
2, 2, 1, 2, 1, 2,
2, 1, 2, 1, 1, 1,
1, 2, 1, 2, 2, 2,
2, 2, 1, 1, 2, 3
};
// clang-format on
for (vtkm::Id i = 0; i < component.GetNumberOfValues(); ++i)
{
VTKM_TEST_ASSERT(component.ReadPortal().Get(i) == componentExpected[size_t(i)],
"Components has unexpected value.");
}
}
};

39
vtkm/worklet/testing/UnitTestPointerJumping.cxx Normal file
View File

@ -0,0 +1,39 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#include <vtkm/cont/ArrayHandleCounting.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/worklet/connectivities/UnionFind.h>
void TestLinear()
{
const vtkm::Id N = 100;
auto counting = vtkm::cont::make_ArrayHandleCounting(-1, 1, N - 1);
vtkm::cont::ArrayHandle<vtkm::Id> parents;
vtkm::cont::ArrayCopy(counting, parents);
parents.WritePortal().Set(0, 0);
vtkm::cont::Invoker invoker;
invoker(vtkm::worklet::connectivity::PointerJumping{}, parents);
VTKM_TEST_ASSERT(vtkm::cont::testing::test_equal_ArrayHandles(
vtkm::cont::ArrayHandleConstant<vtkm::Id>(0, N - 1), parents));
}
void TestPointerJumping()
{
TestLinear();
}
int UnitTestPointerJumping(int argc, char* argv[])
{
return vtkm::cont::testing::Testing::Run(TestPointerJumping, argc, argv);
}