zfp 1d worklet.

This commit is contained in:
Mark Kim 2018-12-07 16:02:16 -05:00
parent ece5913869
commit 8a9bfcba1b
14 changed files with 371 additions and 50 deletions

@ -43,6 +43,8 @@ public:
// 1D uniform datasets.
vtkm::cont::DataSet Make1DUniformDataSet0();
vtkm::cont::DataSet Make1DUniformDataSet1();
vtkm::cont::DataSet Make1DUniformDataSet2();
// 1D explicit datasets.
vtkm::cont::DataSet Make1DExplicitDataSet0();
@ -113,6 +115,34 @@ inline vtkm::cont::DataSet MakeTestDataSet::Make1DUniformDataSet1()
return dataSet;
}
//Make a simple 1D, 16 cell uniform dataset.
inline vtkm::cont::DataSet MakeTestDataSet::Make1DUniformDataSet2()
{
vtkm::cont::DataSetBuilderUniform dsb;
vtkm::Id dims = 16;
vtkm::cont::DataSet dataSet = dsb.Create(dims);
vtkm::cont::DataSetFieldAdd dsf;
const vtkm::Id nVerts = 256;
vtkm::Float64 pointvar[dims];
vtkm::Float64 dx = vtkm::Float64(4.0 * vtkm::Pi()) / vtkm::Float64(dims - 1);
vtkm::Id idx = 0;
for (vtkm::Id x = 0; x < dims; ++x)
{
vtkm::Float64 cx = vtkm::Float64(x) * dx - 2.0 * vtkm::Pi();
vtkm::Float64 cv = vtkm::Sin(cx);
pointvar[idx] = cv;
idx++;
}
dsf.AddPointField(dataSet, "pointvar", pointvar, nVerts);
return dataSet;
}
inline vtkm::cont::DataSet MakeTestDataSet::Make1DExplicitDataSet0()
{
const int nVerts = 5;

@ -22,8 +22,8 @@
#define vtk_m_filter_ZFPDecompressor1D_h
#include <vtkm/filter/FilterField.h>
#include <vtkm/worklet/ZFPCompressor.h>
#include <vtkm/worklet/ZFPDecompress.h>
#include <vtkm/worklet/ZFP1DDecompress.h>
#include <vtkm/worklet/zfp/ZFPTools.h>
namespace vtkm
{
@ -66,7 +66,7 @@ public:
private:
vtkm::Float64 rate;
vtkm::worklet::ZFPDecompressor decompressor;
vtkm::worklet::ZFP1DDecompressor decompressor;
};
template <>

@ -75,10 +75,9 @@ inline VTKM_CONT vtkm::cont::DataSet ZFPDecompressor1D::DoExecute(
hasCellFields = true;
}
}
const vtkm::Id3 dim(field.GetNumberOfValues(), 1, 1);
vtkm::cont::ArrayHandle<vtkm::Float32, StorageType> decompress;
decompressor.Decompress(field, decompress, rate, dim);
vtkm::cont::ArrayHandle<vtkm::Float64, StorageType> decompress;
decompressor.Decompress(field, decompress, rate, field.GetNumberOfValues());
vtkm::cont::DataSet dataset;
vtkm::cont::Field decompressField(

@ -36,29 +36,43 @@
namespace vtkm_ut_zfp_filter
{
void TestZFP1DFilter()
void TestZFP1DFilter(vtkm::Float64 rate)
{
vtkm::cont::testing::MakeTestDataSet testDataSet;
vtkm::cont::DataSet dataset = testDataSet.Make1DUniformDataSet0();
vtkm::cont::DataSet dataset = testDataSet.Make1DUniformDataSet2();
auto dynField = dataset.GetField("pointvar").GetData();
vtkm::cont::ArrayHandle<vtkm::Float64> field =
dynField.Cast<vtkm::cont::ArrayHandle<vtkm::Float64>>();
auto oport = field.GetPortalControl();
vtkm::filter::ZFPCompressor1D compressor;
vtkm::filter::ZFPDecompressor1D decompressor;
compressor.SetActiveField("pointvar");
compressor.SetRate(4);
compressor.SetRate(rate);
auto compressed = compressor.Execute(dataset);
decompressor.SetActiveField("compressed");
decompressor.SetRate(4);
decompressor.SetRate(rate);
auto decompress = decompressor.Execute(compressed);
dynField = decompress.GetField("decompressed").GetData();
;
field = dynField.Cast<vtkm::cont::ArrayHandle<vtkm::Float64>>();
auto port = field.GetPortalControl();
for (int i = 0; i < field.GetNumberOfValues(); i++)
{
std::cout << oport.Get(i) << " " << port.Get(i) << " " << oport.Get(i) - port.Get(i)
<< std::endl;
;
}
}
void TestZFP2DFilter()
void TestZFP2DFilter(vtkm::Float64 rate)
{
@ -75,13 +89,13 @@ void TestZFP2DFilter()
vtkm::filter::ZFPDecompressor2D decompressor;
compressor.SetActiveField("pointvar");
compressor.SetRate(4);
compressor.SetRate(rate);
auto compressed = compressor.Execute(dataset);
decompressor.SetActiveField("compressed");
decompressor.SetRate(4);
decompressor.SetRate(rate);
auto decompress = decompressor.Execute(compressed);
dynField = decompress.GetField("decompressed").GetData();
;
@ -98,8 +112,8 @@ void TestZFP2DFilter()
void TestZFPFilter()
{
//TestZFP1DFilter();
TestZFP2DFilter();
TestZFP1DFilter(4);
//TestZFP2DFilter(4);
}
} // anonymous namespace

@ -87,6 +87,7 @@ set(headers
ZFPCompressor.h
ZFPDecompress.h
ZFP1DCompressor.h
ZFP1DDecompress.h
ZFP2DCompressor.h
ZFP2DDecompress.h
)

@ -0,0 +1,126 @@
//============================================================================
// 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.
//
// Copyright 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-NA0003525 with NTESS,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_worklet_zfp_1d_decompressor_h
#define vtk_m_worklet_zfp_1d_decompressor_h
#include <vtkm/Math.h>
#include <vtkm/cont/Algorithm.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/ArrayHandleConstant.h>
#include <vtkm/cont/ArrayHandleCounting.h>
#include <vtkm/cont/AtomicArray.h>
#include <vtkm/cont/Timer.h>
#include <vtkm/cont/testing/MakeTestDataSet.h>
#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/worklet/zfp/ZFPDecode1.h>
#include <vtkm/worklet/zfp/ZFPTools.h>
using ZFPWord = vtkm::UInt64;
#include <stdio.h>
namespace vtkm
{
namespace worklet
{
namespace detail
{
} // namespace detail
class ZFP1DDecompressor
{
public:
template <typename Scalar>
void Decompress(const vtkm::cont::ArrayHandle<vtkm::Int64>& encodedData,
vtkm::cont::ArrayHandle<Scalar>& output,
const vtkm::Float64 requestedRate,
vtkm::Id dims)
{
//DataDumpb(data, "uncompressed");
zfp::ZFPStream stream;
constexpr vtkm::Int32 topoDims = 1;
;
vtkm::Float64 actualRate = stream.SetRate(requestedRate, topoDims, vtkm::Float64());
std::cout << "ArraySize " << encodedData.GetNumberOfValues() << "\n";
std::cout << "Array dims " << dims << "\n";
std::cout << "requested rate " << requestedRate << " actual rate " << actualRate << "\n";
std::cout << "MinBits " << stream.minbits << "\n";
// Check to see if we need to increase the block sizes
// in the case where dim[x] is not a multiple of 4
vtkm::Id paddedDims = dims;
// ensure that we have block sizes
// that are a multiple of 4
if (paddedDims % 4 != 0)
paddedDims += 4 - dims % 4;
constexpr vtkm::Id four = 4;
vtkm::Id totalBlocks = (paddedDims / four);
std::cout << "Padded dims " << paddedDims << "\n";
size_t outbits = detail::CalcMem1d(paddedDims, stream.minbits);
std::cout << "Total output bits " << outbits << "\n";
vtkm::Id outsize = outbits / sizeof(ZFPWord);
std::cout << "Output size " << outsize << "\n";
output.Allocate(dims);
// hopefully this inits/allocates the mem only on the device
//
//vtkm::cont::ArrayHandleConstant<vtkm::Int64> zero(0, outsize);
//vtkm::cont::Algorithm::Copy(zero, output);
//
using Timer = vtkm::cont::Timer<vtkm::cont::DeviceAdapterTagSerial>;
{
Timer timer;
vtkm::cont::ArrayHandleCounting<vtkm::Id> one(0, 1, 1);
vtkm::worklet::DispatcherMapField<detail::MemTransfer> dis;
dis.Invoke(one, output);
dis.Invoke(one, encodedData);
vtkm::Float64 time = timer.GetElapsedTime();
std::cout << "Copy scalars " << time << "\n";
}
// launch 1 thread per zfp block
vtkm::cont::ArrayHandleCounting<vtkm::Id> blockCounter(0, 1, totalBlocks);
Timer timer;
vtkm::worklet::DispatcherMapField<zfp::Decode1> decompressDispatcher(
zfp::Decode1(dims, paddedDims, stream.maxbits));
decompressDispatcher.Invoke(blockCounter, output, encodedData);
vtkm::Float64 time = timer.GetElapsedTime();
size_t total_bytes = output.GetNumberOfValues() * sizeof(vtkm::Float64);
vtkm::Float64 gB = vtkm::Float64(total_bytes) / (1024. * 1024. * 1024.);
vtkm::Float64 rate = gB / time;
std::cout << "Decompress time " << time << " sec\n";
std::cout << "Decompress rate " << rate << " GB / sec\n";
DataDump(output, "decompressed");
}
};
} // namespace worklet
} // namespace vtkm
#endif // vtk_m_worklet_zfp_1d_decompressor_h

@ -120,9 +120,9 @@ public:
vtkm::Float64 rate = gB / time;
std::cout << "Decompress time " << time << " sec\n";
std::cout << "Decompress rate " << rate << " GB / sec\n";
DataDumpb(output, "decompressed");
DataDump(output, "decompressed");
}
};
} // namespace worklet
} // namespace vtkm
#endif // vtk_m_worklet_zfp_2d_compressor_h
#endif // vtk_m_worklet_zfp_2d_decompressor_h

@ -77,33 +77,6 @@ namespace detail
} // namespace detail
template <typename T>
T* GetVTKMPointerb(vtkm::cont::ArrayHandle<T>& handle)
{
typedef typename vtkm::cont::ArrayHandle<T> HandleType;
typedef typename HandleType::template ExecutionTypes<vtkm::cont::DeviceAdapterTagSerial>::Portal
PortalType;
typedef typename vtkm::cont::ArrayPortalToIterators<PortalType>::IteratorType IteratorType;
IteratorType iter =
vtkm::cont::ArrayPortalToIterators<PortalType>(handle.GetPortalControl()).GetBegin();
return &(*iter);
}
template <typename T>
void DataDumpb(vtkm::cont::ArrayHandle<T> handle, std::string fileName)
{
T* ptr = GetVTKMPointerb(handle);
vtkm::Id osize = handle.GetNumberOfValues();
FILE* fp = fopen(fileName.c_str(), "wb");
;
if (fp != NULL)
{
fwrite(ptr, sizeof(T), osize, fp);
}
fclose(fp);
}
class ZFPDecompressor
{
@ -180,7 +153,7 @@ public:
vtkm::Float64 rate = gB / time;
std::cout << "Decompress time " << time << " sec\n";
std::cout << "Decompress rate " << rate << " GB / sec\n";
DataDumpb(output, "decompressed");
DataDump(output, "decompressed");
}
};
} // namespace worklet

@ -21,6 +21,9 @@
#include <vtkm/worklet/ZFPCompressor.h>
#include <vtkm/worklet/ZFPDecompress.h>
#include <vtkm/worklet/ZFP1DCompressor.h>
#include <vtkm/worklet/ZFP1DDecompress.h>
#include <vtkm/worklet/ZFP2DCompressor.h>
#include <vtkm/worklet/ZFP2DDecompress.h>
@ -30,7 +33,42 @@
#include <iostream>
using Handle64 = vtkm::cont::ArrayHandle<vtkm::Float64>;
template <typename Scalar>
void Test1D(int rate)
{
std::cout << "Testing ZFP 1d:" << std::endl;
vtkm::Id dims = 16;
vtkm::cont::testing::MakeTestDataSet testDataSet;
vtkm::cont::DataSet dataset = testDataSet.Make1DUniformDataSet2();
auto dynField = dataset.GetField("pointvar").GetData();
vtkm::worklet::ZFP1DCompressor compressor;
vtkm::worklet::ZFP1DDecompressor decompressor;
if (dynField.IsSameType(Handle64()))
{
Handle64 field = dynField.Cast<Handle64>();
vtkm::cont::ArrayHandle<Scalar> handle;
const vtkm::Id size = field.GetNumberOfValues();
handle.Allocate(size);
auto fPortal = field.GetPortalControl();
auto hPortal = handle.GetPortalControl();
for (vtkm::Id i = 0; i < size; ++i)
{
hPortal.Set(i, static_cast<Scalar>(fPortal.Get(i)));
}
auto compressed = compressor.Compress(handle, rate, dims);
vtkm::cont::ArrayHandle<Scalar> decoded;
decompressor.Decompress(compressed, decoded, rate, dims);
auto port = decoded.GetPortalControl();
for (int i = 0; i < decoded.GetNumberOfValues(); i++)
{
std::cout << port.Get(i) << std::endl;
}
}
}
template <typename Scalar>
void Test2D(int rate)
{
@ -102,8 +140,9 @@ void Test3D(int rate)
void TestZFP()
{
Test3D<vtkm::Float64>(4);
Test2D<vtkm::Float64>(4);
// Test3D<vtkm::Float64>(4);
// Test2D<vtkm::Float64>(4);
Test1D<vtkm::Float64>(4);
//Test3D<vtkm::Float32>(4);
//Test3D<vtkm::Int64>(4);
//Test3D<vtkm::Int32>(4);

@ -22,6 +22,7 @@ set(headers
ZFPBlockWriter.h
ZFPCodec.h
ZFPDecode.h
ZFPDecode1.h
ZFPDecode2.h
ZFPDecode3.h
ZFPEncode.h

@ -0,0 +1,122 @@
#ifndef vtk_m_worklet_zfp_decode1_h
#define vtk_m_worklet_zfp_decode1_h
#include <vtkm/Types.h>
#include <vtkm/internal/ExportMacros.h>
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/zfp/ZFPBlockWriter.h>
#include <vtkm/worklet/zfp/ZFPDecode.h>
#include <vtkm/worklet/zfp/ZFPFunctions.h>
#include <vtkm/worklet/zfp/ZFPStructs.h>
#include <vtkm/worklet/zfp/ZFPTypeInfo.h>
namespace vtkm
{
namespace worklet
{
namespace zfp
{
template <typename Scalar, typename PortalType>
VTKM_EXEC inline void ScatterPartial1(const Scalar* q,
PortalType& scalars,
const vtkm::Id dims,
vtkm::Id offset,
vtkm::Int32 nx)
{
vtkm::Id x, y;
for (x = 0; x < nx; x++, offset++, q++)
{
scalars.Set(offset, *q);
}
}
template <typename Scalar, typename PortalType>
VTKM_EXEC inline void Scatter1(const Scalar* q,
PortalType& scalars,
const vtkm::Id dims,
vtkm::Id offset)
{
for (vtkm::Id x = 0; x < 4; x++, ++offset)
{
scalars.Set(offset, *q++);
} // x
}
struct Decode1 : public vtkm::worklet::WorkletMapField
{
protected:
vtkm::Id Dims; // field dims
vtkm::Id PaddedDims; // dims padded to a multiple of zfp block size
vtkm::Id ZFPDims; // zfp block dims
vtkm::UInt32 MaxBits; // bits per zfp block
public:
Decode1(const vtkm::Id dims, const vtkm::Id paddedDims, const vtkm::UInt32 maxbits)
: Dims(dims)
, PaddedDims(paddedDims)
, MaxBits(maxbits)
{
ZFPDims = PaddedDims / 4;
}
using ControlSignature = void(FieldIn<>, WholeArrayOut<>, WholeArrayIn<> bitstream);
using ExecutionSignature = void(_1, _2, _3);
template <typename InputScalarPortal, typename BitstreamPortal>
VTKM_EXEC void operator()(const vtkm::Id blockIdx,
InputScalarPortal& scalars,
BitstreamPortal& stream) const
{
using Scalar = typename InputScalarPortal::ValueType;
constexpr vtkm::Int32 BlockSize = 4;
Scalar fblock[BlockSize];
// clear
for (vtkm::Int32 i = 0; i < BlockSize; ++i)
{
fblock[i] = static_cast<Scalar>(0);
}
zfp::zfp_decode<BlockSize>(fblock, MaxBits, blockIdx, stream);
//for(int i = 0; i < BlockSize; ++i)
//{
// std::cout<<" "<<fblock[i];
//}
std::cout << "\n";
vtkm::Id zfpBlock;
zfpBlock = blockIdx % ZFPDims;
vtkm::Id logicalStart = zfpBlock * vtkm::Id(4);
//std::cout<<"Block ID "<<blockIdx<<"\n";
//std::cout<<"ZFP Block "<<zfpBlock<<"\n";
//std::cout<<"logicalStart Start "<<logicalStart<<"\n";
// get the offset into the field
//vtkm::Id offset = (zfpBlock[2]*4*ZFPDims[1] + zfpBlock[1] * 4)*ZFPDims[0] * 4 + zfpBlock[0] * 4;
//std::cout<<"ZFP block offset "<<offset<<"\n";
bool partial = false;
if (logicalStart + 4 > Dims)
partial = true;
//std::cout<<"Dims "<<Dims<<"\n";
if (partial)
{
const vtkm::Int32 nx =
logicalStart + 4 > Dims ? vtkm::Int32(Dims - logicalStart) : vtkm::Int32(4);
//std::cout<<"Partial block "<<logicalStart<<" offset "<<offset<<"\n";
//std::cout<<"Nx "<<nx<<" "<<ny<<" "<<nz<<"\n";
ScatterPartial1(fblock, scalars, Dims, logicalStart, nx);
}
else
{
//std::cout<<"FULL block "<<zfpBlock<<"\n";
Scatter1(fblock, scalars, Dims, logicalStart);
}
}
};
}
}
} // namespace vtkm::worklet::zfp
#endif

@ -214,6 +214,21 @@ inline VTKM_EXEC void fwd_xform<vtkm::Int32, 16>(vtkm::Int32* p)
for (x = 0; x < 4; x++)
fwd_lift<vtkm::Int32, 4>(p + 1 * x);
}
template <>
inline VTKM_EXEC void fwd_xform<vtkm::Int64, 4>(vtkm::Int64* p)
{
/* transform along x */
fwd_lift<vtkm::Int64, 1>(p);
}
template <>
inline VTKM_EXEC void fwd_xform<vtkm::Int32, 4>(vtkm::Int32* p)
{
/* transform along x */
fwd_lift<vtkm::Int32, 1>(p);
}
template <vtkm::Int32 BlockSize, typename PortalType, typename Int>
VTKM_EXEC void encode_block(BlockWriter<BlockSize, PortalType>& stream,
vtkm::Int32 maxbits,

@ -28,7 +28,7 @@ VTKM_EXEC inline void GatherPartial1(Scalar* q,
int sx)
{
vtkm::Id x;
for (x = 0; x < nx; x++, offset += 1)
for (x = 0; x < nx; x++, offset += sx)
q[x] = scalars.Get(offset);
PadBlock(q, nx, 1);
}

@ -81,10 +81,10 @@ size_t CalcMem2d(const vtkm::Id2 dims, const int bits_per_block)
return alloc_size * sizeof(ZFPWord);
}
size_t CalcMem1d(const vtkm::Id2 dims, const int bits_per_block)
size_t CalcMem1d(const vtkm::Id dims, const int bits_per_block)
{
constexpr size_t vals_per_block = 4;
const size_t size = dims[0];
const size_t size = dims;
size_t total_blocks = size / vals_per_block;
constexpr size_t bits_per_word = sizeof(ZFPWord) * 8;
const size_t total_bits = bits_per_block * total_blocks;
@ -123,6 +123,7 @@ void DataDump(vtkm::cont::ArrayHandle<T> handle, std::string fileName)
}
} // namespace worklet
} // namespace vtkm
#endif // vtk_m_worklet_zfp_tools_h