1#ifndef GRID_HPP
2#define GRID_HPP
3
4
5#include "config.h"
6#include "util/cuda_launch.hpp"
7#include <boost/shared_array.hpp>
8#include <vector>
9#include <initializer_list>
10#include <array>
11#include "memory/memory.hpp"
12#include "Space/Shape/Box.hpp"
13#include "grid_key.hpp"
14#include <iostream>
15#include "util/mathutil.hpp"
16#include "iterators/stencil_type.hpp"
17
18
19// Box need the definition of grid_key_dx_r
20#define HARDWARE 1
21
22
23struct No_check
24{
25 template<typename SparseGrid_type>
26 __device__ __host__ bool check(SparseGrid_type & sggt,unsigned int dataBlockPos, unsigned int offset)
27 {
28 return true;
29 }
30};
31
32template<unsigned int dim, typename T>
33struct Box_check
34{
35 Box<dim,T> box;
36
37 template<typename T2>
38 explicit Box_check(Box<dim,T2> & box)
39 :box(box)
40 {}
41
42 template<typename SparseGridGpu_type>
43 __device__ __host__ bool check(SparseGridGpu_type & sggt,unsigned int dataBlockId, unsigned int offset)
44 {
45 auto key = sggt.getCoord(dataBlockId,offset);
46
47 return box.isInsideKey(key);
48 }
49};
50
51
52/*! \brief Class to check if the edge can be created or not
53 *
54 * Class to check if the edge can be created or not, in this case no check is implemented
55 *
56 */
57
58class NoCheck
59{
60public:
61 /*! \brief No check is performed
62 *
63 * No check is performed
64 *
65 * \param v_id Vertex id
66 * \param sz Size limit for the vertex id
67 *
68 */
69 static bool valid(size_t v_id, size_t sz)
70 {
71 return true;
72 }
73};
74
75/*! \brief Class to check if the edge can be created or not
76 *
77 * Class to check if the edge can be created or not, in this case no check is implemented
78 *
79 */
80
81class CheckExistence
82{
83public:
84 /*! \brief Check if vertex exist
85 *
86 * Check if exist
87 *
88 * \param v_id Vertex id
89 * \param sz Size limit for the vertex id
90 *
91 * \return true if exist
92 *
93 */
94 static bool valid(size_t v_id, size_t sz)
95 {
96 return v_id < sz;
97 }
98};
99
100template<unsigned int dim>
101struct ite_gpu
102{
103#ifdef CUDA_GPU
104
105 dim3 thr;
106 dim3 wthr;
107
108 grid_key_dx<dim,int> start;
109 grid_key_dx<dim,int> stop;
110
111 size_t nblocks()
112 {
113 return wthr.x * wthr.y * wthr.z;
114 }
115
116 size_t nthrs()
117 {
118 return thr.x * thr.y * thr.z;
119 }
120
121#endif
122};
123
124//! Declaration grid_sm
125template<unsigned int N, typename T> class grid_sm;
126
127template<unsigned int dim, typename T2, typename T>
128ite_gpu<dim> getGPUIterator_impl(const grid_sm<dim,T2> & g1, const grid_key_dx<dim,T> & key1, const grid_key_dx<dim,T> & key2, size_t n_thr = 1024);
129
130//! Declaration print_warning_on_adjustment
131template <unsigned int dim, typename linearizer> class print_warning_on_adjustment;
132
133//! Declaration grid_key_dx_iterator_sub
134template<unsigned int dim,typename stencil=no_stencil, typename linearizer = grid_sm<dim,void>, typename warn=print_warning_on_adjustment<dim,linearizer>> class grid_key_dx_iterator_sub;
135
136/*! \brief class that store the information of the grid like number of point on each direction and
137 * define the index linearization by stride
138 *
139 * \param N dimensionality
140 * \param T type of object is going to store the grid
141 *
142 */
143template<unsigned int N, typename T>
144class grid_sm
145{
146 //! Box enclosing the grid
147 Box<N,size_t> box;
148
149 //! total number of the elements in the grid
150 size_t size_tot;
151
152 //! size of the grid
153 size_t sz[N];
154
155 //! size of the grid on each stride (used for linearization)
156 size_t sz_s[N];
157
158 /*! \brief Initialize the basic structure
159 *
160 * Initialize the basic structure
161 *
162 * \param sz vector that store the size of the grid on each
163 * dimensions
164 *
165 */
166
167 inline void Initialize(const size_t sz)
168 {
169 //! Initialize the basic structure for each dimension
170 sz_s[0] = sz;
171 this->sz[0] = sz;
172
173 // set the box
174 box.setHigh(0,sz);
175 box.setLow(0,0);
176
177 for (size_t i = 1 ; i < N ; i++)
178 {
179 /* coverity[dead_error_begin] */
180 sz_s[i] = sz*sz_s[i-1];
181 this->sz[i] = sz;
182
183 // set the box
184 box.setHigh(i,sz);
185 box.setLow(i,0);
186 }
187 }
188
189 /*! \brief Initialize the basic structure
190 *
191 * Initialize the basic structure
192 *
193 * \param sz vector that store the size of the grid on each
194 * dimensions
195 *
196 */
197
198 inline void Initialize(const size_t (& sz)[N])
199 {
200 //! Initialize the basic structure for each dimension
201 sz_s[0] = sz[0];
202 this->sz[0] = sz[0];
203
204 // set the box
205 box.setHigh(0,sz[0]);
206 box.setLow(0,0);
207
208 for (size_t i = 1 ; i < N ; i++)
209 {
210 /* coverity[dead_error_begin] */
211 sz_s[i] = sz[i]*sz_s[i-1];
212 this->sz[i] = sz[i];
213
214 // set the box
215 box.setHigh(i,sz[i]);
216 box.setLow(i,0);
217 }
218 }
219
220 /*! \brief Initialize the basic structure
221 *
222 * Produce a grid of size 0 on each dimension
223 *
224 */
225
226 inline void Initialize()
227 {
228 //! Initialize the basic structure for each dimension
229 sz_s[0] = 0;
230 this->sz[0] = 0;
231
232 // set the box
233 box.setHigh(0,0);
234 box.setLow(0,0);
235
236 for (size_t i = 1 ; i < N ; i++)
237 {
238 /* coverity[dead_error_begin] */
239 sz_s[i] = sz[i]*sz_s[i-1];
240
241 // set the box
242 box.setHigh(i,sz[i]);
243 box.setLow(i,0);
244 }
245 }
246
247 /*! \brief linearize an arbitrary set of index
248 *
249 * linearize an arbitrary set of index
250 *
251 */
252 template<typename a, typename ...lT>
253 __device__ __host__ inline mem_id Lin_impl(a v,lT...t) const
254 {
255 return v*sz_s[sizeof...(t)-1] + Lin_impl(t...);
256 }
257
258 //! Linearize a set of index
259 template<typename a> __device__ __host__ inline mem_id Lin_impl(a v) const
260 {
261 return v;
262 }
263
264public:
265
266
267 /*! \brief Return the box enclosing the grid
268 *
269 * \return the box
270 *
271 */
272 inline Box<N,size_t> getBox() const
273 {
274 return box;
275 }
276
277 /*! \brief Return the box enclosing the grid
278 *
279 * While getBox return as P2 the size of the grid
280 * getBoxKey return the size - 1 equivalent to the maximum valid point
281 * that does not overflow the grid
282 *
283 * \return the box
284 *
285 */
286 inline const Box<N,size_t> getBoxKey() const
287 {
288 Box<N,size_t> bx;
289
290 for (size_t i = 0 ; i < N ; i++)
291 {
292 bx.setLow(i,box.getLow(i));
293 bx.setHigh(i,box.getHigh(i) - 1);
294 }
295
296 return bx;
297 }
298
299 /*! \brief Reset the dimension of the grid
300 *
301 * \param dims store on each dimension the size of the grid
302 *
303 */
304 inline void setDimensions(const size_t (& dims)[N])
305 {
306 Initialize(dims);
307 size_tot = totalSize(dims);
308 }
309
310 /*! \brief Default constructor
311 *
312 * It produce a grid of size 0 on each dimension
313 *
314 */
315
316 inline grid_sm()
317 :size_tot(0)
318 {
319 // Initialize sz
320 for (size_t i = 0 ; i < N ; i++)
321 {sz[i] = 0;}
322
323 Initialize();
324 }
325
326 /*! \brief construct a grid from another grid
327 *
328 * \param g grid info
329 *
330 * construct a grid from another grid, type can be different
331 *
332 */
333
334 template<typename S> inline grid_sm(const grid_sm<N,S> & g)
335 {
336 size_tot = g.size_tot;
337
338 for (size_t i = 0 ; i < N ; i++)
339 {
340 sz[i] = g.sz[i];
341 sz_s[i] = g.sz_s[i];
342 }
343 }
344
345 // Static element to calculate total size
346
347 inline size_t totalSize(const size_t sz)
348 {
349 size_t tSz = 1;
350
351 for (size_t i = 0 ; i < N ; i++)
352 {
353 tSz *= sz;
354 }
355
356 return tSz;
357 }
358
359 // Static element to calculate total size
360
361 inline size_t totalSize(const size_t (& sz)[N])
362 {
363 size_t tSz = 1;
364
365 for (size_t i = 0 ; i < N ; i++)
366 {
367 tSz *= sz[i];
368 }
369
370 return tSz;
371 }
372
373
374 /*! \brief Construct a grid of a specified size
375 *
376 * Construct a grid of a specified size
377 *
378 * \param sz is an array that contain the size of the grid on each dimension
379 *
380 */
381
382 inline grid_sm(const size_t & sz)
383 : size_tot(totalSize(sz))
384 {
385 Initialize(sz);
386 }
387
388 /*! \brief Construct a grid of a specified size
389 *
390 * Construct a grid of a specified size
391 *
392 * \param sz is an array that contain the size of the grid on each dimension
393 *
394 */
395
396 inline grid_sm(const size_t (& sz)[N])
397 : size_tot(totalSize(sz))
398 {
399 Initialize(sz);
400 }
401
402 /*! \brief Linearization of the grid_key_dx with a specified shift
403 *
404 * \tparam check class that check the linearization, if this check fail the function return -1
405 * \param gk grid_key_dx to linearize
406 * \param sum_id shift on each dimension
407 *
408 * \return The linearization of the gk key shifted by c, or -1 if the check fail
409 */
410
411 template<typename check=NoCheck, typename ids_type>
412 inline mem_id LinId(const grid_key_dx<N,ids_type> & gk, const char sum_id[N]) const
413 {
414 mem_id lid;
415
416 // Check the sum produce a valid key
417
418 if (check::valid(gk.get(0) + sum_id[0],sz[0]) == false)
419 return -1;
420
421 lid = gk.get(0) + sum_id[0];
422
423
424 for (mem_id i = 1 ; i < N ; i++)
425 {
426 // Check the sum produce a valid key
427
428 if (check::valid(gk.get(i) + sum_id[i],sz[i]) == false)
429 return -1;
430
431 lid += (gk.get(i) + sum_id[i]) * sz_s[i-1];
432 }
433
434 return lid;
435 }
436
437 /*! \brief Linearization of the grid_key_dx with a specified shift
438 *
439 * \tparam check class that check the linearization, if this check fail the function return -1
440 * \param gk grid_key_dx to linearize
441 * \param sum_id shift on each dimension
442 * \param bc boundary conditions
443 *
444 * \return The linearization of the gk key shifted by c, or -1 if the check fail
445 */
446
447 template<typename check=NoCheck,typename ids_type>
448 inline mem_id LinId(const grid_key_dx<N,ids_type> & gk, const char sum_id[N], const size_t (&bc)[N]) const
449 {
450 mem_id lid;
451
452 // Check the sum produce a valid key
453
454 if (bc[0] == NON_PERIODIC)
455 {
456 if (check::valid(gk.get(0) + sum_id[0],sz[0]) == false)
457 return -1;
458
459 lid = gk.get(0) + sum_id[0];
460 }
461 else
462 {
463 lid = openfpm::math::positive_modulo(gk.get(0) + sum_id[0],sz[0]);
464 }
465
466 for (mem_id i = 1 ; i < N ; i++)
467 {
468 // Check the sum produce a valid key
469
470 /* coverity[dead_error_line] */
471 if (bc[i] == NON_PERIODIC)
472 {
473 if (check::valid(gk.get(i) + sum_id[i],sz[i]) == false)
474 return -1;
475
476 lid += (gk.get(i) + sum_id[i]) * sz_s[i-1];
477 }
478 else
479 {
480 lid += (openfpm::math::positive_modulo(gk.get(i) + sum_id[i],sz[i])) * sz_s[i-1];
481 }
482 }
483
484 return lid;
485 }
486
487 /*! \brief Linearization of the set of indexes
488 *
489 * Linearization of the set of indexes, it spit out a number that is just the 1D linearization.
490 * In this case is the linearization of N index
491 *
492 * \param k set of indexes to linearize
493 *
494 */
495 inline mem_id LinIdPtr(size_t * k) const
496 {
497 mem_id lid = k[0];
498 for (mem_id i = 1 ; i < N ; i++)
499 {
500 lid += k[i] * sz_s[i-1];
501 }
502
503 return lid;
504 }
505
506 /*! \brief Linearization of the grid_key_dx
507 *
508 * Linearization of the grid_key_dx given a key, it spit out a number that is just the 1D linearization
509 * of the key. In this case is the linearization of N index
510 *
511 * \param k grid key to access the element on the grid
512 *
513 */
514
515 __device__ __host__ inline mem_id LinId(const size_t (& k)[N]) const
516 {
517 mem_id lid = k[0];
518 for (mem_id i = 1 ; i < N ; i++)
519 {
520 /* coverity[dead_error_line] */
521 lid += k[i] * sz_s[i-1];
522 }
523
524 return lid;
525 }
526
527 /*! \brief Linearization of the grid_key_dx
528 *
529 * Linearization of the grid_key_dx given a key, it spit out a number that is just the 1D linearization
530 * of the key. In this case is the linearization of N index
531 *
532 * \param gk grid key to access the element of the grid
533 *
534 */
535
536 template<typename ids_type> __device__ __host__ inline mem_id LinId(const grid_key_dx<N,ids_type> & gk) const
537 {
538 mem_id lid = gk.get(0);
539 for (mem_id i = 1 ; i < N ; i++)
540 {
541 /* coverity[dead_error_line] */
542 lid += gk.get(i) * sz_s[i-1];
543 }
544
545 return lid;
546 }
547
548 /*! \brief linearize an arbitrary set of index
549 *
550 * linearize an arbitrary set of index
551 *
552 */
553 template<typename a, typename ...lT, typename enabler = typename std::enable_if<sizeof...(lT) == N-1>::type >
554 __device__ __host__ inline mem_id Lin(a v,lT...t) const
555 {
556 return v*sz_s[sizeof...(t)-1] + Lin_impl(t...);
557 }
558
559 //! Construct
560
561 /*! \brief inversion of the linearization of the grid_key_dx
562 *
563 * \param id of the object
564 * \return key of the grid that id identify
565 *
566 */
567 __device__ __host__ inline grid_key_dx<N> InvLinId(mem_id id) const
568 {
569 // Inversion of linearize
570
571 grid_key_dx<N> gk;
572
573 for (mem_id i = 0 ; i < N ; i++)
574 {
575 gk.set_d(i,id % sz[i]);
576 id /= sz[i];
577 }
578
579 return gk;
580 }
581
582
583 /*! \brief Linearization of an array of mem_id (long int)
584 *
585 * Linearization of an array of mem_id, it spit out a number that is just the 1D linearization
586 * of the key. In this case is the linearization of N index
587 *
588 * \param id an array of mem_id index
589 *
590 */
591
592 //#pragma openfpm layout(get)
593 inline mem_id LinId(mem_id * id) const
594 {
595 mem_id lid = 0;
596 lid += id[0];
597 for (mem_id i = 1 ; i < N ; i++)
598 {
599 lid += id[i] * sz_s[i-1];
600 }
601
602 return lid;
603 }
604
605 //! Destructor
606 ~grid_sm() {};
607
608 /*! \brief Return the size of the grid
609 *
610 * Return the size of the grid
611 *
612 * \return the size of the grid
613 *
614 */
615 __device__ __host__ inline size_t size() const
616 {
617 return size_tot;
618 };
619
620 /*! \brief Copy the grid from another grid
621 *
622 * \param g grid from witch to copy
623 *
624 */
625
626 __device__ __host__ inline grid_sm<N,T> & operator=(const grid_sm<N,T> & g)
627 {
628 size_tot = g.size_tot;
629
630 for (size_t i = 0 ; i < N ; i++)
631 {
632 sz[i] = g.sz[i];
633 sz_s[i] = g.sz_s[i];
634 }
635
636 box = g.box;
637
638 return *this;
639 }
640
641 /*! \brief Check if the two grid_sm are the same
642 *
643 * \param g element to check
644 *
645 * \return true if they are the same
646 *
647 */
648
649 inline bool operator==(const grid_sm<N,T> & g)
650 {
651 for (size_t i = 0 ; i < N ; i++)
652 {
653 if (sz[i] != g.sz[i])
654 return false;
655 }
656
657#ifdef SE_CLASS1
658
659 if (size_tot != g.size_tot)
660 return false;
661
662 for (size_t i = 0 ; i < N ; i++)
663 {
664 if (sz_s[i] != g.sz_s[i])
665 return false;
666 }
667
668#endif
669 return true;
670 }
671
672 /*! \brief Check if the two grid_sm are the same
673 *
674 * \param g element to check
675 *
676 */
677
678 inline bool operator!=(const grid_sm<N,T> & g)
679 {
680 return ! this->operator==(g);
681 }
682
683 /**
684 *
685 * Get the stride-size of the grid on the direction i
686 *
687 * [Example] on a grid 16*16*16 it return 16,256
688 *
689 * \param i direction
690 * \return the size on the direction i
691 *
692 */
693
694 inline size_t size_s(unsigned int i) const
695 {
696 return sz_s[i];
697 }
698
699 /**
700 *
701 * Get the size of the grid on the direction i
702 *
703 * \param i direction
704 * \return the size on the direction i
705 *
706 */
707
708 __device__ __host__ inline size_t size(unsigned int i) const
709 {
710 return sz[i];
711 }
712
713 /*! \brief Return the size of the grid as an array
714 *
715 * \return get the size of the grid as an array
716 *
717 */
718 inline const size_t (& getSize() const)[N]
719 {
720 return sz;
721 }
722
723 /*! \brief Return a sub-grid iterator
724 *
725 * Return a sub-grid iterator, to iterate through the grid
726 *
727 * \param start start point
728 * \param stop stop point
729 *
730 */
731 inline grid_key_dx_iterator_sub<N> getSubIterator(const grid_key_dx<N> & start, const grid_key_dx<N> & stop) const
732 {
733 return grid_key_dx_iterator_sub<N>(*this,start,stop);
734 }
735
736#ifdef CUDA_GPU
737
738 /*! \brief Get an iterator for the GPU
739 *
740 * \param start starting point
741 * \param stop end point
742 *
743 */
744 template<typename T2>
745 struct ite_gpu<N> getGPUIterator(const grid_key_dx<N,T2> & key1, const grid_key_dx<N,T2> & key2, size_t n_thr = 1024) const
746 {
747 return getGPUIterator_impl<N>(*this,key1,key2,n_thr);
748 }
749
750 /*! \brief Get an iterator for the GPU
751 *
752 * \param start starting point
753 * \param stop end point
754 *
755 */
756 struct ite_gpu<N> getGPUIterator(size_t n_thr = 1024) const
757 {
758 grid_key_dx<N> k1;
759 grid_key_dx<N> k2;
760
761 for (size_t i = 0 ; i < N ; i++)
762 {
763 k1.set_d(i,0);
764 k2.set_d(i,size(i));
765 }
766
767 return getGPUIterator_impl<N>(*this,k1,k2,n_thr);
768 }
769
770#endif
771
772 /*! \brief swap the grid_sm informations
773 *
774 * \param g grid to swap
775 *
776 */
777 inline void swap(grid_sm<N,T> & g)
778 {
779 size_t tmp = size_tot;
780 size_tot = g.size_tot;
781 g.size_tot = tmp;
782
783 for (size_t i = 0 ; i < N ; i++)
784 {
785 tmp = sz[i];
786 sz[i] = g.sz[i];
787 g.sz[i] = tmp;
788
789 tmp = sz_s[i];
790 sz_s[i] = g.sz_s[i];
791 g.sz_s[i] = tmp;
792 }
793 }
794
795 /*! \brief Produce a string from the object
796 *
797 * \return string
798 *
799 */
800 std::string toString() const
801 {
802 std::stringstream str;
803
804 for (size_t i = 0 ; i < N ; i++)
805 str << "sz[" << i << "]=" << size(i) << " ";
806
807 return str.str();
808 }
809
810 //! It simply mean that all the classes grid are friend of all its specialization
811 template <unsigned int,typename> friend class grid_sm;
812};
813
814
815template<unsigned int dim, typename T2, typename T>
816ite_gpu<dim> getGPUIterator_impl(const grid_sm<dim,T2> & g1, const grid_key_dx<dim,T> & key1, const grid_key_dx<dim,T> & key2, const size_t n_thr)
817{
818 size_t tot_work = 1;
819 for (size_t i = 0 ; i < dim ; i++)
820 {tot_work *= key2.get(i) - key1.get(i) + 1;}
821
822 size_t n = (tot_work <= n_thr)?openfpm::math::round_big_2(tot_work):n_thr;
823
824 // Work to do
825 ite_gpu<dim> ig;
826
827 if (tot_work == 0)
828 {
829 ig.thr.x = 0;
830 ig.thr.y = 0;
831 ig.thr.z = 0;
832
833 ig.wthr.x = 0;
834 ig.wthr.y = 0;
835 ig.wthr.z = 0;
836
837 return ig;
838 }
839
840 ig.thr.x = 1;
841 ig.thr.y = 1;
842 ig.thr.z = 1;
843
844 int dir = 0;
845
846 while (n != 1)
847 {
848 if (dir % 3 == 0)
849 {ig.thr.x = ig.thr.x << 1;}
850 else if (dir % 3 == 1)
851 {ig.thr.y = ig.thr.y << 1;}
852 else if (dir % 3 == 2)
853 {ig.thr.z = ig.thr.z << 1;}
854
855 n = n >> 1;
856
857 dir++;
858 dir %= dim;
859 }
860
861 if (dim >= 1)
862 {ig.wthr.x = (key2.get(0) - key1.get(0) + 1) / ig.thr.x + (((key2.get(0) - key1.get(0) + 1)%ig.thr.x != 0)?1:0);}
863
864 if (dim >= 2)
865 {ig.wthr.y = (key2.get(1) - key1.get(1) + 1) / ig.thr.y + (((key2.get(1) - key1.get(1) + 1)%ig.thr.y != 0)?1:0);}
866 else
867 {ig.wthr.y = 1;}
868
869 if (dim >= 3)
870 {
871 // Roll the other dimensions on z
872 ig.wthr.z = 1;
873 for (size_t i = 2 ; i < dim ; i++)
874 {ig.wthr.z *= (key2.get(i) - key1.get(i) + 1) / ig.thr.z + (((key2.get(i) - key1.get(i) + 1)%ig.thr.z != 0)?1:0);}
875 }
876 else
877 {ig.wthr.z = 1;}
878
879 // crop if wthr == 1
880
881 if (dim >= 1 && ig.wthr.x == 1)
882 {ig.thr.x = (key2.get(0) - key1.get(0) + 1);}
883
884 if (dim >= 2 && ig.wthr.y == 1)
885 {ig.thr.y = key2.get(1) - key1.get(1) + 1;}
886
887 if (dim == 3 && ig.wthr.z == 1)
888 {ig.thr.z = key2.get(2) - key1.get(2) + 1;}
889
890 for (size_t i = 0 ; i < dim ; i++)
891 {
892 ig.start.set_d(i,key1.get(i));
893 ig.stop.set_d(i,key2.get(i));
894 }
895
896 return ig;
897}
898
899
900/*! \brief Emulate grid_key_dx with runtime dimensionality
901 *
902 * Emulate grid_key_dx with runtime dimensionality
903 *
904 */
905
906class grid_key_dx_r
907{
908 size_t dim;
909
910public:
911
912 /*! \brief Get the dimensionality of the key
913 *
914 * Get the dimensionality of the key
915 *
916 */
917
918 size_t getDim()
919 {
920 return dim;
921 }
922
923 /*! \brief constructor from another key
924 *
925 * \param key
926 *
927 */
928 grid_key_dx_r(grid_key_dx_r & key)
929 :dim(key.dim)
930 {
931 // Allocate the key
932 k = new mem_id[dim];
933
934 // Copy the key
935 for(unsigned int i = 0 ; i < dim ; i++)
936 {
937 k[i] = key.k[i];
938 }
939 }
940
941 /*! \brief constructor
942 *
943 * constructor
944 *
945 * \param dim Dimensionality
946 *
947 */
948 grid_key_dx_r(size_t dim)
949 :dim(dim)
950 {
951 // Allocate the key
952 k = new mem_id[dim];
953 }
954
955 ~grid_key_dx_r()
956 {
957 delete [] k;
958 }
959
960 //! set the grid key from a list of numbers
961 template<typename a, typename ...T>void set(a v, T...t)
962 {
963 k[dim-1] = v;
964 invert_assign(t...);
965 }
966
967 /*! \brief get the i index
968 *
969 * Get the i index
970 *
971 * \param i index to get
972 *
973 * \return the index value
974 *
975 */
976 mem_id get(size_t i)
977 {
978 return k[i];
979 }
980
981 /*! \brief Set the i index
982 *
983 * Set the i index
984 *
985 * \param i index to set
986 * \param id value to set
987 *
988 */
989 void set_d(size_t i, mem_id id)
990 {
991 k[i] = id;
992 }
993
994 //! structure that store all the index
995 mem_id * k;
996
997private:
998
999 /*! \brief Recursively invert the assignment
1000 *
1001 * Recursively invert the assignment at compile-time
1002 *
1003 */
1004 template<typename a, typename ...T>void invert_assign(a v,T...t)
1005 {
1006 k[sizeof...(T)] = v;
1007 invert_assign(t...);
1008 }
1009
1010 template<typename a, typename ...T>void invert_assign(a v)
1011 {
1012 k[0] = v;
1013 }
1014
1015 void invert_assign()
1016 {
1017 }
1018
1019};
1020
1021
1022/**
1023 *
1024 * Iterate through the elements (i1,i2,....,in) with i1 ... in unsigned integers
1025 * with the following constrain (i1>i2>......>in)
1026 *
1027 */
1028
1029class Iterator_g_const
1030{
1031 size_t dim;
1032
1033 //! size of the grid (the grid is assumed a square so equal on each dimension)
1034 size_t sz;
1035
1036 // Actual grid_key position
1037 grid_key_dx_r gk;
1038
1039public:
1040
1041 /*! \brief Get the dimensionality of the iterator
1042 *
1043 * Get the dimensionality of the iterator
1044 *
1045 */
1046
1047 size_t getDim()
1048 {
1049 return dim;
1050 }
1051
1052 /*! \brief Constructor
1053 *
1054 * \param n Dimensionality (how many i1 ... in you have)
1055 * \param sz Size of the grid on all dimensions range of the value i1 ... in can assume
1056 *
1057 */
1058
1059 Iterator_g_const(size_t n, size_t sz)
1060 :dim(n),sz(sz),gk(n)
1061 {
1062 // fill gk with the first grid element that satisfied the constrain: 0,1,2,3... dim
1063
1064 for (size_t i = 0 ; i < dim ; i++)
1065 {
1066 gk.set_d(i,dim-i-1);
1067 }
1068 }
1069
1070 /*! \brief Get the next element
1071 *
1072 * Get the next element
1073 *
1074 * \return the next grid_key
1075 *
1076 */
1077
1078 Iterator_g_const & operator++()
1079 {
1080 //! increment the first index
1081
1082 gk.set_d(0,gk.get(0)+1);
1083
1084 //! check the overflow of all the index with exception of the last dimensionality
1085
1086 unsigned int i = 0;
1087 for ( ; i < dim-1 ; i++)
1088 {
1089 size_t id = gk.get(i);
1090 if (id >= sz)
1091 {
1092 // ! overflow, increment the next index
1093
1094 id = gk.get(i+1);
1095 if (id+i+2 >= sz)
1096 {
1097 // there is no-way to produce a valid key
1098 // there is not need to check the previous index
1099 // overflow i+1
1100
1101 gk.set_d(i+1,sz);
1102 }
1103 else
1104 {
1105
1106 // reinitialize the previous index
1107
1108 for (unsigned int s = 0 ; s <= i+1 ; s++)
1109 {
1110 gk.set_d(i+1-s,id+1+s);
1111 }
1112 }
1113 }
1114 else
1115 {
1116 break;
1117 }
1118 }
1119
1120 return *this;
1121 }
1122
1123 /*! \brief Check if there is the next element
1124 *
1125 * Check if there is the next element
1126 *
1127 * \return true if there is the next, false otherwise
1128 *
1129 */
1130
1131 bool isNext()
1132 {
1133 // If dimensionless return immediately
1134 if (dim == 0)
1135 return false;
1136
1137 if (gk.get(dim-1) < static_cast<mem_id>(sz-dim+1))
1138 {
1139 //! we did not reach the end of the grid
1140
1141 return true;
1142 }
1143
1144 //! we reach the end of the grid
1145 return false;
1146 }
1147
1148 /*! \brief Return the actual key
1149 *
1150 * Return the actual key
1151 *
1152 * \return The actual key that identify with the set of index
1153 *
1154 */
1155
1156 grid_key_dx_r & get()
1157 {
1158 return gk;
1159 }
1160};
1161
1162#endif
1163