fml  0.1-0
Fused Matrix Library
arraytools.hpp
1 // This file is part of arraytools which is released under the Boost Software
2 // License, Version 1.0. See accompanying file LICENSE or copy at
3 // https://www.boost.org/LICENSE_1_0.txt
4 
5 #ifndef SPVEC_ARRAYTOOLS_H
6 #define SPVEC_ARRAYTOOLS_H
7 #pragma once
8 
9 
10 #include <algorithm>
11 #include <cfloat>
12 #include <cmath>
13 #include <cstdlib>
14 #include <cstring>
15 #include <stdexcept>
16 
17 #ifdef __CUDACC__
18 #include <cublas.h>
19 #endif
20 
21 
22 namespace arraytools
23 {
30  template <typename T>
31  static inline void alloc(const size_t len, T **x)
32  {
33  *x = (T*) std::malloc(len*sizeof(T));
34  }
35 
36 
37 
44  template <typename T>
45  static inline void zero_alloc(const size_t len, T **x)
46  {
47  *x = (T*) std::calloc(len, sizeof(T));
48  }
49 
50 
51 
59  template <typename T>
60  static inline void realloc(const size_t len, T **x)
61  {
62  void *realloc_ptr = std::realloc(*x, len*sizeof(T));
63  if (realloc_ptr == NULL)
64  std::free(*x);
65 
66  *x = (T*) realloc_ptr;
67  }
68 
69 
70 
76  template <typename T>
77  static inline void free(T *x)
78  {
79  if (x)
80  std::free(x);
81  }
82 
83 
84 
93  template <typename T>
94  static inline void copy(const size_t len, const T *src, T *dst)
95  {
96  std::memcpy(dst, src, len*sizeof(*src));
97  }
98 
99  template <typename SRC, typename DST>
100  static inline void copy(const size_t len, const SRC *src, DST *dst)
101  {
102  #pragma omp for simd
103  for (size_t i=0; i<len; i++)
104  dst[i] = (DST) src[i];
105  }
106 
107 
108 
115  template <typename T>
116  static inline void zero(const size_t len, T *x)
117  {
118  std::memset(x, 0, len*sizeof(*x));
119  }
120 
121 
122 
123  namespace
124  {
125  static const int ALLOC_FAILED = -1;
126  static const int ALLOC_OK = 0;
127 
128 
129 
130  static inline int check_alloc()
131  {
132  return ALLOC_OK;
133  }
134 
135  template <typename T>
136  static inline int check_alloc(T *x)
137  {
138  if (x == NULL)
139  return ALLOC_FAILED;
140  else
141  return ALLOC_OK;
142  }
143 
144  template <typename T, typename... VAT>
145  static inline int check_alloc(T *x, VAT... vax)
146  {
147  return check_alloc(x) + check_alloc(vax ...);
148  }
149 
150 
151 
152  static inline void vafree(){}
153 
154  template <typename T, typename... VAT>
155  static inline void vafree(T *x, VAT... vax)
156  {
157  arraytools::free(x);
158  vafree(vax ...);
159  }
160  }
161 
162 
163 
171  template <typename T, typename... VAT>
172  static inline void check_allocs(T *x, VAT... vax)
173  {
174  int check = check_alloc(x) + check_alloc(vax ...);
175 
176  if (check != ALLOC_OK)
177  {
178  arraytools::free(x);
179  vafree(vax ...);
180 
181  throw std::bad_alloc();
182  }
183  }
184 
185 
186 
187  namespace fltcmp
188  {
189  namespace
190  {
191  // eps ~ sqrt(precision)
192  #ifdef __CUDACC__
193  const __half eps_hlf = 1e-2; // machine epsilon 0x1p-11 ~ 2.2e-2
194  #endif
195  const float eps_flt = 1e-4;
196  const double eps_dbl = 1e-8;
197 
198  // smallest representible number with 1 as the leading bit of the significand
199  #ifdef __CUDACC__
200  const __half min_normal_hlf = 0x1p-14; // (2^(5-1)-1)-1 = 14
201  #endif
202  const float min_normal_flt = 0x1p-126;
203  const double min_normal_dbl = 0x1p-1022;
204 
205  #ifdef __CUDACC__
206  const __half max_hlf = 1>>15;
207  #endif
208  const float max_flt = FLT_MAX;
209  const double max_dbl = DBL_MAX;
210 
211 
212 
213  #ifdef __CUDACC__
214  static inline __half fabsh(const __half x)
215  {
216  return __float2half(fabsf(__half2float(x)));
217  }
218  #endif
219 
220 
221 
222  static inline bool subnormal(const float x)
223  {
224  return (x < min_normal_flt);
225  }
226 
227  static inline bool subnormal(const double x)
228  {
229  return (x < min_normal_dbl);
230  }
231 
232 
233 
234  template <typename T>
235  static inline bool abseq(const T x, const T y)
236  {
237  float fx = (float) x;
238  float fy = (float) y;
239  return fabsf(fx - fy) < (eps_flt * min_normal_flt);
240  }
241 
242  static inline bool abseq(const float x, const float y)
243  {
244  return fabsf(x - y) < (eps_flt * min_normal_flt);
245  }
246 
247  static inline bool abseq(const double x, const double y)
248  {
249  return fabs(x - y) < (eps_dbl * min_normal_dbl);
250  }
251 
252 
253 
254  template <typename T>
255  static inline bool releq(const T x, const T y)
256  {
257  float fx = (float) x;
258  float fy = (float) y;
259  return fabsf(fx - fy) / std::min(fabsf(fx)+fabsf(fy), max_flt) < eps_flt;
260  }
261 
262  static inline bool releq(const float x, const float y)
263  {
264  return fabsf(x - y) / std::min(fabsf(x)+fabsf(y), max_flt) < eps_flt;
265  }
266 
267  static inline bool releq(const double x, const double y)
268  {
269  return fabs(x - y) / std::min(fabs(x)+fabs(y), max_dbl) < eps_dbl;
270  }
271  }
272 
273 
274 
275  // modified from https://floating-point-gui.de/errors/comparison/
276  template <typename REAL>
277  static inline bool eq(const REAL x, const REAL y)
278  {
279  if (x == y)
280  return true;
281  else if (x == 0.f || y == 0.f || subnormal(fabs(x)+fabs(y)))
282  return abseq(x, y);
283  else
284  return releq(x, y);
285  }
286 
287 
288 
289  static inline bool eq(const int x, const int y)
290  {
291  return (x == y);
292  }
293 
294 
295 
296  template <typename REAL>
297  static inline bool eq(const REAL x, const int y)
298  {
299  return eq(x, (REAL) y);
300  }
301 
302  template <typename REAL>
303  static inline bool eq(const int x, const REAL y)
304  {
305  return eq((REAL) x, y);
306  }
307 
308 
309 
310  static inline bool eq(const float x, double y)
311  {
312  return eq(x, (float) y);
313  }
314 
315  static inline bool eq(const double x, const float y)
316  {
317  return eq((float) x, y);
318  }
319  }
320 
321 
322 
331  template <typename TA, typename TB>
332  static inline size_t cmp_firstmiss(const size_t len, const TA *a, const TB *b)
333  {
334  size_t ret = len;
335 
336  #pragma omp simd reduction(min:ret)
337  for (size_t i=0; i<len; i++)
338  {
339  if (!fltcmp::eq(a[i], b[i]))
340  ret = i;
341  }
342 
343  return ret;
344  }
345 
346 
347 
356  template <typename TA, typename TB>
357  static inline bool cmp(const size_t len, const TA *a, const TB *b)
358  {
359  size_t fm = cmp_firstmiss(len, a, b);
360  return (fm == len ? true : false);
361  }
362 }
363 
364 
365 #endif
Definition: arraytools.hpp:22