VTK-m  2.0
ZFPEncode.h
Go to the documentation of this file.
1 //============================================================================
2 // Copyright (c) Kitware, Inc.
3 // All rights reserved.
4 // See LICENSE.txt for details.
5 //
6 // This software is distributed WITHOUT ANY WARRANTY; without even
7 // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
8 // PURPOSE. See the above copyright notice for more information.
9 //============================================================================
10 #ifndef vtk_m_worklet_zfp_encode_h
11 #define vtk_m_worklet_zfp_encode_h
12 
13 #include <vtkm/Types.h>
18 
19 namespace vtkm
20 {
21 namespace worklet
22 {
23 namespace zfp
24 {
25 
26 template <typename Scalar>
28 {
29  switch (n)
30  {
31  case 0:
32  p[0 * s] = 0;
33  /* FALLTHROUGH */
34  case 1:
35  p[1 * s] = p[0 * s];
36  /* FALLTHROUGH */
37  case 2:
38  p[2 * s] = p[1 * s];
39  /* FALLTHROUGH */
40  case 3:
41  p[3 * s] = p[0 * s];
42  /* FALLTHROUGH */
43  default:
44  break;
45  }
46 }
47 
48 template <vtkm::Int32 N, typename FloatType>
50 {
51  FloatType maxVal = 0;
52  for (vtkm::Int32 i = 0; i < N; ++i)
53  {
54  maxVal = vtkm::Max(maxVal, vtkm::Abs(vals[i]));
55  }
56 
57  if (maxVal > 0)
58  {
59  vtkm::Int32 exponent;
60  vtkm::Frexp(maxVal, &exponent);
61  /* clamp exponent in case x is denormal */
62  return vtkm::Max(exponent, 1 - get_ebias<FloatType>());
63  }
64  return -get_ebias<FloatType>();
65 }
66 
67 // maximum number of bit planes to encode
69 {
70  return vtkm::Min(maxprec, vtkm::Max(0, maxexp - minexp + 8));
71 }
72 
73 template <typename Scalar>
74 inline VTKM_EXEC Scalar quantize(Scalar x, vtkm::Int32 e)
75 {
76  return vtkm::Ldexp(x, (CHAR_BIT * (vtkm::Int32)sizeof(Scalar) - 2) - e);
77 }
78 
79 template <typename Int, typename Scalar, vtkm::Int32 BlockSize>
80 inline VTKM_EXEC void fwd_cast(Int* iblock, const Scalar* fblock, vtkm::Int32 emax)
81 {
82  Scalar s = quantize<Scalar>(1, emax);
83  for (vtkm::Int32 i = 0; i < BlockSize; ++i)
84  {
85  iblock[i] = static_cast<Int>(s * fblock[i]);
86  }
87 }
88 
89 template <typename Int, vtkm::Int32 S>
90 inline VTKM_EXEC void fwd_lift(Int* p)
91 {
92  Int x, y, z, w;
93  x = *p;
94  p += S;
95  y = *p;
96  p += S;
97  z = *p;
98  p += S;
99  w = *p;
100  p += S;
101 
102  /*
103  ** non-orthogonal transform
104  ** ( 4 4 4 4) (x)
105  ** 1/16 * ( 5 1 -1 -5) (y)
106  ** (-4 4 4 -4) (z)
107  ** (-2 6 -6 2) (w)
108  */
109  x += w;
110  x >>= 1;
111  w -= x;
112  z += y;
113  z >>= 1;
114  y -= z;
115  x += z;
116  x >>= 1;
117  z -= x;
118  w += y;
119  w >>= 1;
120  y -= w;
121  w += y >> 1;
122  y -= w >> 1;
123 
124  p -= S;
125  *p = w;
126  p -= S;
127  *p = z;
128  p -= S;
129  *p = y;
130  p -= S;
131  *p = x;
132 }
133 
134 template <typename Int, typename UInt>
135 inline VTKM_EXEC UInt int2uint(const Int x);
136 
137 template <>
138 inline VTKM_EXEC vtkm::UInt64 int2uint<vtkm::Int64, vtkm::UInt64>(const vtkm::Int64 x)
139 {
140  return (static_cast<vtkm::UInt64>(x) + (vtkm::UInt64)0xaaaaaaaaaaaaaaaaull) ^
141  (vtkm::UInt64)0xaaaaaaaaaaaaaaaaull;
142 }
143 
144 template <>
145 inline VTKM_EXEC vtkm::UInt32 int2uint<vtkm::Int32, vtkm::UInt32>(const vtkm::Int32 x)
146 {
147  return (static_cast<vtkm::UInt32>(x) + (vtkm::UInt32)0xaaaaaaaau) ^ (vtkm::UInt32)0xaaaaaaaau;
148 }
149 
150 
151 
152 template <typename UInt, typename Int, vtkm::Int32 BlockSize>
153 inline VTKM_EXEC void fwd_order(UInt* ublock, const Int* iblock)
154 {
155  const zfp::ZFPCodec<BlockSize> codec;
156  for (vtkm::Int32 i = 0; i < BlockSize; ++i)
157  {
158  vtkm::UInt8 idx = codec.CodecLookup(i);
159  ublock[i] = int2uint<Int, UInt>(iblock[idx]);
160  }
161 }
162 
163 template <typename Int, vtkm::Int32 BlockSize>
164 inline VTKM_EXEC void fwd_xform(Int* p);
165 
166 template <>
167 inline VTKM_EXEC void fwd_xform<vtkm::Int64, 64>(vtkm::Int64* p)
168 {
169  vtkm::UInt32 x, y, z;
170  /* transform along x */
171  for (z = 0; z < 4; z++)
172  for (y = 0; y < 4; y++)
173  fwd_lift<vtkm::Int64, 1>(p + 4 * y + 16 * z);
174  /* transform along y */
175  for (x = 0; x < 4; x++)
176  for (z = 0; z < 4; z++)
177  fwd_lift<vtkm::Int64, 4>(p + 16 * z + 1 * x);
178  /* transform along z */
179  for (y = 0; y < 4; y++)
180  for (x = 0; x < 4; x++)
181  fwd_lift<vtkm::Int64, 16>(p + 1 * x + 4 * y);
182 }
183 
184 template <>
185 inline VTKM_EXEC void fwd_xform<vtkm::Int32, 64>(vtkm::Int32* p)
186 {
187  vtkm::UInt32 x, y, z;
188  /* transform along x */
189  for (z = 0; z < 4; z++)
190  for (y = 0; y < 4; y++)
191  fwd_lift<vtkm::Int32, 1>(p + 4 * y + 16 * z);
192  /* transform along y */
193  for (x = 0; x < 4; x++)
194  for (z = 0; z < 4; z++)
195  fwd_lift<vtkm::Int32, 4>(p + 16 * z + 1 * x);
196  /* transform along z */
197  for (y = 0; y < 4; y++)
198  for (x = 0; x < 4; x++)
199  fwd_lift<vtkm::Int32, 16>(p + 1 * x + 4 * y);
200 }
201 
202 template <>
203 inline VTKM_EXEC void fwd_xform<vtkm::Int64, 16>(vtkm::Int64* p)
204 {
205  vtkm::UInt32 x, y;
206  /* transform along x */
207  for (y = 0; y < 4; y++)
208  fwd_lift<vtkm::Int64, 1>(p + 4 * y);
209  /* transform along y */
210  for (x = 0; x < 4; x++)
211  fwd_lift<vtkm::Int64, 4>(p + 1 * x);
212 }
213 
214 template <>
215 inline VTKM_EXEC void fwd_xform<vtkm::Int32, 16>(vtkm::Int32* p)
216 {
217  vtkm::UInt32 x, y;
218  /* transform along x */
219  for (y = 0; y < 4; y++)
220  fwd_lift<vtkm::Int32, 1>(p + 4 * y);
221  /* transform along y */
222  for (x = 0; x < 4; x++)
223  fwd_lift<vtkm::Int32, 4>(p + 1 * x);
224 }
225 
226 template <>
227 inline VTKM_EXEC void fwd_xform<vtkm::Int64, 4>(vtkm::Int64* p)
228 {
229  /* transform along x */
230  fwd_lift<vtkm::Int64, 1>(p);
231 }
232 
233 template <>
234 inline VTKM_EXEC void fwd_xform<vtkm::Int32, 4>(vtkm::Int32* p)
235 {
236  /* transform along x */
237  fwd_lift<vtkm::Int32, 1>(p);
238 }
239 
240 template <vtkm::Int32 BlockSize, typename PortalType, typename Int>
242  vtkm::Int32 maxbits,
243  vtkm::Int32 maxprec,
244  Int* iblock)
245 {
246  using UInt = typename zfp_traits<Int>::UInt;
247 
248  fwd_xform<Int, BlockSize>(iblock);
249 
250  UInt ublock[BlockSize];
251  fwd_order<UInt, Int, BlockSize>(ublock, iblock);
252 
253  vtkm::UInt32 intprec = CHAR_BIT * (vtkm::UInt32)sizeof(UInt);
254  vtkm::UInt32 kmin =
255  intprec > (vtkm::UInt32)maxprec ? intprec - static_cast<vtkm::UInt32>(maxprec) : 0;
256  vtkm::UInt32 bits = static_cast<vtkm::UInt32>(maxbits);
257  vtkm::UInt32 i, m;
258  vtkm::UInt32 n = 0;
259  vtkm::UInt64 x;
260  /* encode one bit plane at a time from MSB to LSB */
261  for (vtkm::UInt32 k = intprec; bits && k-- > kmin;)
262  {
263  /* step 1: extract bit plane #k to x */
264  x = 0;
265  for (i = 0; i < BlockSize; i++)
266  {
267  x += (vtkm::UInt64)((ublock[i] >> k) & 1u) << i;
268  }
269  /* step 2: encode first n bits of bit plane */
270  m = vtkm::Min(n, bits);
271  bits -= m;
272  x = stream.write_bits(x, m);
273  /* step 3: unary run-length encode remainder of bit plane */
274  for (; n < BlockSize && bits && (bits--, stream.write_bit(!!x)); x >>= 1, n++)
275  {
276  for (; n < BlockSize - 1 && bits && (bits--, !stream.write_bit(x & 1u)); x >>= 1, n++)
277  {
278  }
279  }
280  }
281 }
282 
283 
284 template <vtkm::Int32 BlockSize, typename Scalar, typename PortalType>
285 inline VTKM_EXEC void zfp_encodef(Scalar* fblock,
286  vtkm::Int32 maxbits,
287  vtkm::UInt32 blockIdx,
288  PortalType& stream)
289 {
290  using Int = typename zfp::zfp_traits<Scalar>::Int;
291  zfp::BlockWriter<BlockSize, PortalType> blockWriter(stream, maxbits, vtkm::Id(blockIdx));
292  vtkm::Int32 emax = zfp::MaxExponent<BlockSize, Scalar>(fblock);
293  // std::cout<<"EMAX "<<emax<<"\n";
294  vtkm::Int32 maxprec =
295  zfp::precision(emax, zfp::get_precision<Scalar>(), zfp::get_min_exp<Scalar>());
296  vtkm::UInt32 e = vtkm::UInt32(maxprec ? emax + zfp::get_ebias<Scalar>() : 0);
297  /* encode block only if biased exponent is nonzero */
298  if (e)
299  {
300 
301  const vtkm::UInt32 ebits = vtkm::UInt32(zfp::get_ebits<Scalar>()) + 1;
302  blockWriter.write_bits(2 * e + 1, ebits);
303 
304  Int iblock[BlockSize];
305  zfp::fwd_cast<Int, Scalar, BlockSize>(iblock, fblock, emax);
306 
307  encode_block<BlockSize>(blockWriter, maxbits - vtkm::Int32(ebits), maxprec, iblock);
308  }
309 }
310 
311 // helpers so we can do partial template instantiation since
312 // the portal type could be on any backend
313 template <vtkm::Int32 BlockSize, typename Scalar, typename PortalType>
315 {
316 };
317 
318 template <vtkm::Int32 BlockSize, typename PortalType>
319 struct ZFPBlockEncoder<BlockSize, vtkm::Float32, PortalType>
320 {
322  vtkm::Int32 maxbits,
323  vtkm::UInt32 blockIdx,
324  PortalType& stream)
325  {
326  zfp_encodef<BlockSize>(fblock, maxbits, blockIdx, stream);
327  }
328 };
329 
330 template <vtkm::Int32 BlockSize, typename PortalType>
331 struct ZFPBlockEncoder<BlockSize, vtkm::Float64, PortalType>
332 {
334  vtkm::Int32 maxbits,
335  vtkm::UInt32 blockIdx,
336  PortalType& stream)
337  {
338  zfp_encodef<BlockSize>(fblock, maxbits, blockIdx, stream);
339  }
340 };
341 
342 template <vtkm::Int32 BlockSize, typename PortalType>
343 struct ZFPBlockEncoder<BlockSize, vtkm::Int32, PortalType>
344 {
346  vtkm::Int32 maxbits,
347  vtkm::UInt32 blockIdx,
348  PortalType& stream)
349  {
350  using Int = typename zfp::zfp_traits<vtkm::Int32>::Int;
351  zfp::BlockWriter<BlockSize, PortalType> blockWriter(stream, maxbits, vtkm::Id(blockIdx));
352  encode_block<BlockSize>(blockWriter, maxbits, get_precision<vtkm::Int32>(), (Int*)fblock);
353  }
354 };
355 
356 template <vtkm::Int32 BlockSize, typename PortalType>
357 struct ZFPBlockEncoder<BlockSize, vtkm::Int64, PortalType>
358 {
359  VTKM_EXEC void encode(vtkm::Int64* fblock,
360  vtkm::Int32 maxbits,
361  vtkm::UInt32 blockIdx,
362  PortalType& stream)
363  {
364  using Int = typename zfp::zfp_traits<vtkm::Int64>::Int;
365  zfp::BlockWriter<BlockSize, PortalType> blockWriter(stream, maxbits, vtkm::Id(blockIdx));
366  encode_block<BlockSize>(blockWriter, maxbits, get_precision<vtkm::Int64>(), (Int*)fblock);
367  }
368 };
369 }
370 }
371 } // namespace vtkm::worklet::zfp
372 #endif
vtkm::worklet::zfp::fwd_order
VTKM_EXEC void fwd_order(UInt *ublock, const Int *iblock)
Definition: ZFPEncode.h:153
vtkm::worklet::zfp::ZFPBlockEncoder
Definition: ZFPEncode.h:314
vtkm::worklet::zfp::ZFPCodec
Definition: ZFPCodec.h:27
VTKM_EXEC
#define VTKM_EXEC
Definition: ExportMacros.h:51
vtkm
Groups connected points that have the same field value.
Definition: Atomic.h:19
Types.h
vtkm::worklet::zfp::quantize
VTKM_EXEC Scalar quantize(Scalar x, vtkm::Int32 e)
Definition: ZFPEncode.h:74
ZFPTypeInfo.h
vtkm::worklet::cellmetrics::FloatType
vtkm::FloatDefault FloatType
Definition: CellAspectFrobeniusMetric.h:50
vtkm::worklet::zfp::ZFPBlockEncoder< BlockSize, vtkm::Int32, PortalType >::encode
VTKM_EXEC void encode(vtkm::Int32 *fblock, vtkm::Int32 maxbits, vtkm::UInt32 blockIdx, PortalType &stream)
Definition: ZFPEncode.h:345
vtkm::Id
vtkm::Int32 Id
Represents an ID (index into arrays).
Definition: Types.h:191
vtkm::worklet::zfp::ZFPBlockEncoder< BlockSize, vtkm::Float64, PortalType >::encode
VTKM_EXEC void encode(vtkm::Float64 *fblock, vtkm::Int32 maxbits, vtkm::UInt32 blockIdx, PortalType &stream)
Definition: ZFPEncode.h:333
vtkm::worklet::zfp::zfp_encodef
VTKM_EXEC void zfp_encodef(Scalar *fblock, vtkm::Int32 maxbits, vtkm::UInt32 blockIdx, PortalType &stream)
Definition: ZFPEncode.h:285
vtkm::worklet::zfp::fwd_xform
VTKM_EXEC void fwd_xform(Int *p)
vtkm::worklet::zfp::fwd_cast
VTKM_EXEC void fwd_cast(Int *iblock, const Scalar *fblock, vtkm::Int32 emax)
Definition: ZFPEncode.h:80
ZFPCodec.h
ExportMacros.h
vtkm::Frexp
VTKM_EXEC_CONT vtkm::Float32 Frexp(vtkm::Float32 x, vtkm::Int32 *exponent)
Decompose floating poing value.
Definition: Math.h:2563
vtkm::worklet::zfp::int2uint
VTKM_EXEC UInt int2uint(const Int x)
vtkm::worklet::zfp::BlockWriter::write_bits
VTKM_EXEC vtkm::UInt64 write_bits(const vtkm::UInt64 &bits, const unsigned int &n_bits)
Definition: ZFPBlockWriter.h:56
vtkm::worklet::zfp::fwd_lift
VTKM_EXEC void fwd_lift(Int *p)
Definition: ZFPEncode.h:90
vtkm::worklet::zfp::ZFPBlockEncoder< BlockSize, vtkm::Int64, PortalType >::encode
VTKM_EXEC void encode(vtkm::Int64 *fblock, vtkm::Int32 maxbits, vtkm::UInt32 blockIdx, PortalType &stream)
Definition: ZFPEncode.h:359
vtkm::worklet::zfp::PadBlock
VTKM_EXEC void PadBlock(Scalar *p, vtkm::UInt32 n, vtkm::UInt32 s)
Definition: ZFPEncode.h:27
vtkm::worklet::zfp::encode_block
VTKM_EXEC void encode_block(BlockWriter< BlockSize, PortalType > &stream, vtkm::Int32 maxbits, vtkm::Int32 maxprec, Int *iblock)
Definition: ZFPEncode.h:241
vtkm::UInt8
uint8_t UInt8
Definition: Types.h:157
vtkm::worklet::zfp::BlockWriter::write_bit
vtkm::UInt32 VTKM_EXEC write_bit(const unsigned int &bit)
Definition: ZFPBlockWriter.h:88
vtkm::worklet::zfp::precision
VTKM_EXEC vtkm::Int32 precision(vtkm::Int32 maxexp, vtkm::Int32 maxprec, vtkm::Int32 minexp)
Definition: ZFPEncode.h:68
vtkm::UInt32
uint32_t UInt32
Definition: Types.h:161
vtkm::worklet::zfp::MaxExponent
VTKM_EXEC vtkm::Int32 MaxExponent(const FloatType *vals)
Definition: ZFPEncode.h:49
vtkm::worklet::zfp::BlockWriter
Definition: ZFPBlockWriter.h:25
vtkm::Float32
float Float32
Definition: Types.h:154
vtkm::Int32
int32_t Int32
Definition: Types.h:160
vtkm::Float64
double Float64
Definition: Types.h:155
vtkm::worklet::zfp::ZFPBlockEncoder< BlockSize, vtkm::Float32, PortalType >::encode
VTKM_EXEC void encode(vtkm::Float32 *fblock, vtkm::Int32 maxbits, vtkm::UInt32 blockIdx, PortalType &stream)
Definition: ZFPEncode.h:321
ZFPBlockWriter.h
vtkm::Ldexp
VTKM_EXEC_CONT vtkm::Float32 Ldexp(vtkm::Float32 x, vtkm::Int32 exponent)
Definition: Math.h:2582
vtkm::worklet::zfp::zfp_traits
Definition: ZFPTypeInfo.h:162