1 #ifndef _MARRAY_UTILITY_HPP_
2 #define _MARRAY_UTILITY_HPP_
3 
4 #include <type_traits>
5 #include <array>
6 #include <vector>
7 #include <utility>
8 #include <iterator>
9 #include <cassert>
10 #include <algorithm>
11 #include <numeric>
12 #include <tuple>
13 #include <cmath>
14 #include <cstddef>
15 #include <string>
16 #include <functional>
17 #include <ostream>
18 
19 #ifdef MARRAY_ENABLE_ASSERTS
20 #define MARRAY_ASSERT(e) assert(e)
21 #else
22 #define MARRAY_ASSERT(e) ((void)0)
23 #endif
24 
25 #include "short_vector.hpp"
26 
27 namespace MArray
28 {
29 
30 /*
31  * The type all_t specifies a range [0,len_i) for an array
32  * dimension i of length len_i (i.e. it selects all of the data along
33  * that dimension).
34  */
all_tMArray::all_t35 struct all_t { constexpr all_t() {} };
bcast_tMArray::bcast_t36 struct bcast_t { constexpr bcast_t() {} };
37 namespace slice
38 {
39     constexpr all_t all;
40     constexpr bcast_t bcast;
41 }
42 
43 #ifndef MARRAY_LEN_TYPE
44 #define MARRAY_LEN_TYPE ptrdiff_t
45 #endif
46 
47 typedef MARRAY_LEN_TYPE len_type;
48 
49 #ifndef MARRAY_STRIDE_TYPE
50 #define MARRAY_STRIDE_TYPE ptrdiff_t
51 #endif
52 
53 typedef MARRAY_STRIDE_TYPE stride_type;
54 
55 #ifndef MARRAY_OPT_NDIM
56 #define MARRAY_OPT_NDIM 6
57 #endif
58 
59 typedef short_vector<len_type,MARRAY_OPT_NDIM> len_vector;
60 typedef short_vector<stride_type,MARRAY_OPT_NDIM> stride_vector;
61 typedef short_vector<unsigned,MARRAY_OPT_NDIM> dim_vector;
62 typedef short_vector<len_type,MARRAY_OPT_NDIM> index_vector;
63 typedef short_vector<unsigned,MARRAY_OPT_NDIM> irrep_vector;
64 template <typename T>
65 using ptr_vector = short_vector<T*,MARRAY_OPT_NDIM>;
66 
67 #ifndef MARRAY_DEFAULT_LAYOUT
68 #define MARRAY_DEFAULT_LAYOUT ROW_MAJOR
69 #endif
70 
71 #define MARRAY_PASTE_(x,y) x##y
72 #define MARRAY_PASTE(x,y) MARRAY_PASTE_(x,y)
73 
74 #define MARRAY_DEFAULT_DPD_LAYOUT_(type) \
75     MARRAY_PASTE(MARRAY_PASTE(type,_),MARRAY_DEFAULT_LAYOUT)
76 
77 #ifndef MARRAY_DEFAULT_DPD_LAYOUT
78 #define MARRAY_DEFAULT_DPD_LAYOUT PREFIX
79 #endif
80 
81 /*
82  * The special value uninitialized is used to construct an array which
83  * does not default- or value-initialize its elements (useful for avoiding
84  * redundant memory operations for scalar types).
85  */
uninitialized_tMArray::uninitialized_t86 struct uninitialized_t { constexpr uninitialized_t() {} };
87 constexpr uninitialized_t uninitialized;
88 
89 /*
90  * Specifies the layout of the array data.
91  */
92 struct layout
93 {
94     int type;
95 
layoutMArray::layout96     constexpr explicit layout(int type) : type(type) {}
97 
operator ==MArray::layout98     bool operator==(layout other) const { return type == other.type; }
operator !=MArray::layout99     bool operator!=(layout other) const { return type != other.type; }
100 };
101 
column_major_layoutMArray::column_major_layout102 struct column_major_layout : layout { constexpr column_major_layout() : layout(0) {} };
103 constexpr column_major_layout COLUMN_MAJOR;
104 
row_major_layoutMArray::row_major_layout105 struct row_major_layout : layout { constexpr row_major_layout() : layout(1) {} };
106 constexpr row_major_layout ROW_MAJOR;
107 
108 constexpr decltype(MARRAY_DEFAULT_LAYOUT) DEFAULT;
109 
110 struct dpd_layout
111 {
112     int type;
113 
dpd_layoutMArray::dpd_layout114     constexpr explicit dpd_layout(int type) : type(type) {}
115 
116     dpd_layout(layout layout);
117 
118     layout base() const;
119 
operator ==MArray::dpd_layout120     bool operator==(dpd_layout other) const { return type == other.type; }
operator !=MArray::dpd_layout121     bool operator!=(dpd_layout other) const { return type != other.type; }
122 };
123 
124 struct balanced_column_major_layout : dpd_layout
balanced_column_major_layoutMArray::balanced_column_major_layout125 { constexpr balanced_column_major_layout() : dpd_layout(0) {} };
126 constexpr balanced_column_major_layout BALANCED_COLUMN_MAJOR;
127 
128 struct balanced_row_major_layout : dpd_layout
balanced_row_major_layoutMArray::balanced_row_major_layout129 { constexpr balanced_row_major_layout() : dpd_layout(1) {} };
130 constexpr balanced_row_major_layout BALANCED_ROW_MAJOR;
131 
132 struct blocked_column_major_layout : dpd_layout
blocked_column_major_layoutMArray::blocked_column_major_layout133 { constexpr blocked_column_major_layout() : dpd_layout(2) {} };
134 constexpr blocked_column_major_layout BLOCKED_COLUMN_MAJOR;
135 
136 struct blocked_row_major_layout : dpd_layout
blocked_row_major_layoutMArray::blocked_row_major_layout137 { constexpr blocked_row_major_layout() : dpd_layout(3) {} };
138 constexpr blocked_row_major_layout BLOCKED_ROW_MAJOR;
139 
140 struct prefix_column_major_layout : dpd_layout
prefix_column_major_layoutMArray::prefix_column_major_layout141 { constexpr prefix_column_major_layout() : dpd_layout(4) {} };
142 constexpr prefix_column_major_layout PREFIX_COLUMN_MAJOR;
143 
144 struct prefix_row_major_layout : dpd_layout
prefix_row_major_layoutMArray::prefix_row_major_layout145 { constexpr prefix_row_major_layout() : dpd_layout(5) {} };
146 constexpr prefix_row_major_layout PREFIX_ROW_MAJOR;
147 
148 constexpr decltype(MARRAY_DEFAULT_DPD_LAYOUT_(BALANCED)) BALANCED;
149 constexpr decltype(MARRAY_DEFAULT_DPD_LAYOUT_(BLOCKED)) BLOCKED;
150 constexpr decltype(MARRAY_DEFAULT_DPD_LAYOUT_(PREFIX)) PREFIX;
151 
dpd_layout(layout layout)152 inline dpd_layout::dpd_layout(layout layout)
153 : type(layout == DEFAULT   ? MARRAY_DEFAULT_DPD_LAYOUT.type :
154        layout == ROW_MAJOR ? MARRAY_PASTE(MARRAY_DEFAULT_DPD_LAYOUT,_ROW_MAJOR).type
155                            : MARRAY_PASTE(MARRAY_DEFAULT_DPD_LAYOUT,_COLUMN_MAJOR).type) {}
156 
base() const157 inline layout dpd_layout::base() const
158 {
159     if (*this == BALANCED_COLUMN_MAJOR ||
160         *this == BLOCKED_COLUMN_MAJOR ||
161         *this == PREFIX_COLUMN_MAJOR)
162     {
163         return COLUMN_MAJOR;
164     }
165     else
166     {
167         return ROW_MAJOR;
168     }
169 }
170 
171 template <typename I> class range_t;
172 
173 namespace detail
174 {
175     template <size_t N>
inverse_permutation(const std::array<unsigned,N> & p)176     std::array<unsigned, N> inverse_permutation(const std::array<unsigned, N>& p)
177     {
178         std::array<unsigned, N> ip;
179         for (unsigned i = 0;i < N;i++) ip[p[i]] = i;
180         return ip;
181     }
182 
inverse_permutation(const dim_vector & p)183     inline dim_vector inverse_permutation(const dim_vector& p)
184     {
185         dim_vector ip(p.size());
186         for (unsigned i = 0;i < p.size();i++) ip[p[i]] = i;
187         return ip;
188     }
189 
190     template <typename...>
191     struct exists {};
192 
193     template <typename T>
194     using decay_t = typename std::decay<T>::type;
195 
196     template <typename T>
197     using remove_cv_t = typename std::remove_cv<T>::type;
198 
199     template <bool Cond, typename T=void>
200     using enable_if_t = typename std::enable_if<Cond,T>::type;
201 
202     template <bool Cond, typename T, typename U>
203     using conditional_t = typename std::conditional<Cond,T,U>::type;
204 
205     template <typename T, typename U=void>
206     using enable_if_integral_t = enable_if_t<std::is_integral<T>::value,U>;
207 
208     template <typename T, typename U=void>
209     using enable_if_not_integral_t = enable_if_t<!std::is_integral<T>::value,U>;
210 
211     template <typename T, typename U=void>
212     using enable_if_const_t = enable_if_t<std::is_const<T>::value,U>;
213 
214     template <typename T, typename U=void>
215     using enable_if_not_const_t = enable_if_t<!std::is_const<T>::value,U>;
216 
217     template <typename T, typename U, typename V=void>
218     using enable_if_assignable_t = enable_if_t<std::is_assignable<T,U>::value,V>;
219 
220     template <typename T, typename U, typename V=void>
221     using enable_if_convertible_t = enable_if_t<std::is_convertible<T,U>::value,V>;
222 
223     template <typename T, typename... Args>
224     struct are_convertible;
225 
226     template <typename T>
227     struct are_convertible<T> : std::true_type {};
228 
229     template <typename T, typename Arg, typename... Args>
230     struct are_convertible<T, Arg, Args...> :
231         conditional_t<std::is_convertible<Arg, T>::value,
232                       are_convertible<T, Args...>,
233                       std::false_type> {};
234 
235     template <typename T, typename... Args>
236     struct are_assignable;
237 
238     template <typename T>
239     struct are_assignable<T> : std::true_type {};
240 
241     template <typename T, typename Arg, typename... Args>
242     struct are_assignable<T, Arg, Args...> :
243         conditional_t<std::is_assignable<T, Arg>::value,
244                       are_assignable<T, Args...>,
245                       std::false_type> {};
246 
247     template <typename T, typename=void>
248     struct is_index_or_slice_helper : std::false_type {};
249 
250     template <typename T>
251     struct is_index_or_slice_helper<T, enable_if_convertible_t<T, int>> : std::true_type {};
252 
253     template <typename I>
254     struct is_index_or_slice_helper<range_t<I>, enable_if_integral_t<I>> : std::true_type {};
255 
256     template <>
257     struct is_index_or_slice_helper<all_t> : std::true_type {};
258 
259     template <>
260     struct is_index_or_slice_helper<bcast_t> : std::true_type {};
261 
262     template <typename T>
263     struct is_index_or_slice : is_index_or_slice_helper<typename std::decay<T>::type> {};
264 
265     template <typename... Args>
266     struct are_indices_or_slices;
267 
268     template<>
269     struct are_indices_or_slices<> : std::true_type {};
270 
271     template <typename Arg, typename... Args>
272     struct are_indices_or_slices<Arg, Args...> :
273         conditional_t<is_index_or_slice<Arg>::value,
274                       are_indices_or_slices<Args...>,
275                       std::false_type> {};
276 
277     template <typename T, typename=void>
278     struct is_container : std::false_type {};
279 
280     template <typename T>
281     struct is_container<T,
282         conditional_t<false,
283                       exists<typename T::value_type,
284                              decltype(std::declval<T>().size()),
285                              decltype(std::declval<T>().begin()),
286                              decltype(std::declval<T>().end())>,
287                       void>>
288     : std::true_type {};
289 
290     template <typename C, typename T, typename=void>
291     struct is_container_of : std::false_type {};
292 
293     template <typename C, typename T>
294     struct is_container_of<C, T, typename std::enable_if<is_container<C>::value>::type>
295     : std::is_assignable<T&, typename C::value_type> {};
296 
297     template <typename C, typename T, typename U=void>
298     using enable_if_container_of_t = typename std::enable_if<is_container_of<C, T>::value,U>::type;
299 
300     template <typename C, typename=void>
301     struct is_container_of_containers : std::false_type {};
302 
303     template <typename C>
304     struct is_container_of_containers<C, enable_if_t<is_container<C>::value>> :
305         is_container<typename C::value_type> {};
306 
307     template <typename C, typename T, typename=void>
308     struct is_container_of_containers_of : std::false_type {};
309 
310     template <typename C, typename T>
311     struct is_container_of_containers_of<C, T, enable_if_t<is_container<C>::value>> :
312         is_container_of<typename C::value_type, T> {};
313 
314     template <typename C, typename T, typename U=void>
315     using enable_if_container_of_containers_of_t =
316         typename std::enable_if<is_container_of_containers_of<C, T>::value,U>::type;
317 
318     template <typename T, typename... Ts>
319     struct are_containers_helper;
320 
321     template <typename T>
322     struct are_containers_helper<T> : is_container<T> {};
323 
324     template <typename T, typename... Ts>
325     struct are_containers_helper
326     : std::conditional<is_container<T>::value,
327                        are_containers_helper<Ts...>,
328                        std::false_type>::type {};
329 
330     template <typename... Ts>
331     struct are_containers;
332 
333     template <>
334     struct are_containers<> : std::true_type {};
335 
336     template <typename... Ts>
337     struct are_containers : are_containers_helper<Ts...> {};
338 
339     template <typename T, typename C, typename... Cs>
340     struct are_containers_of_helper;
341 
342     template <typename T, typename C>
343     struct are_containers_of_helper<T, C> : is_container_of<C, T> {};
344 
345     template <typename T, typename C, typename... Cs>
346     struct are_containers_of_helper
347     : std::conditional<is_container_of<C, T>::value,
348                        are_containers_of_helper<T, Cs...>,
349                        std::false_type>::type {};
350 
351     template <typename T, typename... Cs>
352     struct are_containers_of;
353 
354     template <typename T>
355     struct are_containers_of<T> : std::true_type {};
356 
357     template <typename T, typename... Cs>
358     struct are_containers_of : are_containers_of_helper<T, Cs...> {};
359 
360     template <typename Iterator>
inc_offsets_helper(unsigned i,Iterator)361     void inc_offsets_helper(unsigned i, Iterator) {}
362 
363     template <typename Iterator, typename Offset, typename... Offsets>
inc_offsets_helper(unsigned i,Iterator it,Offset & off0,Offsets &...off)364     void inc_offsets_helper(unsigned i, Iterator it, Offset& off0,
365                             Offsets&... off)
366     {
367         off0 += (*it)[i];
368         inc_offsets_helper(i, ++it, off...);
369     }
370 
371     template <typename Strides, typename... Offsets>
inc_offsets(unsigned i,const Strides & strides,Offsets &...off)372     void inc_offsets(unsigned i, const Strides& strides, Offsets&... off)
373     {
374         inc_offsets_helper(i, strides.begin(), off...);
375     }
376 
377     template <typename Pos, typename Iterator>
dec_offsets_helper(unsigned i,const Pos &,Iterator)378     void dec_offsets_helper(unsigned i, const Pos&, Iterator) {}
379 
380     template <typename Pos, typename Iterator, typename Offset, typename... Offsets>
dec_offsets_helper(unsigned i,const Pos & pos,Iterator it,Offset & off0,Offsets &...off)381     void dec_offsets_helper(unsigned i, const Pos& pos, Iterator it,
382                              Offset& off0, Offsets&... off)
383     {
384         off0 -= pos[i]*(*it)[i];
385         dec_offsets_helper(i, pos, ++it, off...);
386     }
387 
388     template <typename Pos, typename Strides, typename... Offsets>
dec_offsets(unsigned i,const Pos & pos,const Strides & strides,Offsets &...off)389     void dec_offsets(unsigned i, const Pos& pos, const Strides& strides,
390                      Offsets&... off)
391     {
392         dec_offsets_helper(i, pos, strides.begin(), off...);
393     }
394 
395     template <typename Pos, typename Iterator>
move_offsets_helper(const Pos &,Iterator)396     void move_offsets_helper(const Pos&, Iterator) {}
397 
398     template <typename Pos, typename Iterator, typename Offset, typename... Offsets>
move_offsets_helper(const Pos & pos,Iterator it,Offset & off0,Offsets &...off)399     void move_offsets_helper(const Pos& pos, Iterator it,
400                              Offset& off0, Offsets&... off)
401     {
402         for (unsigned i = 0;i < pos.size();i++) off0 += pos[i]*(*it)[i];
403         move_offsets_helper(pos, ++it, off...);
404     }
405 
406     template <typename Pos, typename Strides, typename... Offsets>
move_offsets(const Pos & pos,const Strides & strides,Offsets &...off)407     void move_offsets(const Pos& pos, const Strides& strides,
408                       Offsets&... off)
409     {
410         move_offsets_helper(pos, strides.begin(), off...);
411     }
412 
413     template <typename Iterator>
set_strides_helper(Iterator)414     void set_strides_helper(Iterator) {}
415 
416     template <typename Iterator, typename Stride, typename... Strides>
set_strides_helper(Iterator it,const Stride & stride,const Strides &...strides)417     void set_strides_helper(Iterator it, const Stride& stride, const Strides&... strides)
418     {
419         std::copy(stride.begin(), stride.end(), it->begin());
420         set_strides_helper(++it, strides...);
421     }
422 
423     template <typename Strides_, typename... Strides>
set_strides(Strides_ & strides_,const Strides &...strides)424     void set_strides(Strides_& strides_, const Strides&... strides)
425     {
426         set_strides_helper(strides_.begin(), strides...);
427     }
428 
429     template <typename T, T... S> struct integer_sequence {};
430 
431     template <typename T, typename U, typename V> struct concat_sequences;
432     template <typename T, T... S, T... R>
433     struct concat_sequences<T, integer_sequence<T, S...>, integer_sequence<T, R...>>
434     {
435         typedef integer_sequence<T, S..., (R+sizeof...(S))...> type;
436     };
437 
438     template <typename T, T N, typename=void> struct static_range_helper;
439 
440     template <typename T, T N> struct static_range_helper<T, N, enable_if_t<N==0>>
441     {
442         typedef integer_sequence<T> type;
443     };
444 
445     template <typename T, T N> struct static_range_helper<T, N, enable_if_t<N==1>>
446     {
447         typedef integer_sequence<T, 0> type;
448     };
449 
450     template <typename T, T N> struct static_range_helper<T, N, enable_if_t<(N>1)>>
451     {
452         typedef typename concat_sequences<T, typename static_range_helper<T, (N+1)/2>::type,
453                                              typename static_range_helper<T, N/2>::type>::type type;
454     };
455 
456     template <typename T, T N>
457     using static_range = typename static_range_helper<T, N>::type;
458 }
459 
460 }
461 
462 #endif
463