Intrepid2
Intrepid2_TensorViewIterator.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ************************************************************************
3 //
4 // Intrepid2 Package
5 // Copyright (2007) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Kyungjoo Kim (kyukim@sandia.gov),
38 // Mauro Perego (mperego@sandia.gov), or
39 // Nate Roberts (nvrober@sandia.gov)
40 //
41 // ************************************************************************
42 // @HEADER
43 
49 #ifndef Intrepid2_TensorViewIterator_h
50 #define Intrepid2_TensorViewIterator_h
51 
54 
55 #include "Kokkos_Vector.hpp"
56 #include <vector>
57 
58 namespace Intrepid2
59 {
72  template<class TensorViewType, class ViewType1, class ViewType2 ,typename ScalarType>
74  {
75  public:
76  enum RankCombinationType
77  {
78  DIMENSION_MATCH,
79  TENSOR_PRODUCT,
80  TENSOR_CONTRACTION
81  };
82  using RankCombinationViewType = Kokkos::View<RankCombinationType*, typename TensorViewType::device_type>;
83  protected:
84 
85  ViewIterator<TensorViewType, ScalarType> tensor_view_iterator_;
88 
89  RankCombinationViewType rank_combination_types_;
90  public:
110  KOKKOS_INLINE_FUNCTION
111  TensorViewIterator(TensorViewType tensor_view, ViewType1 view1, ViewType2 view2,
112  RankCombinationViewType rank_combination_types)
113  :
114  tensor_view_iterator_(tensor_view),
115  view1_iterator_(view1),
116  view2_iterator_(view2),
117  rank_combination_types_(rank_combination_types)
118  {
119  // rank_combination_type should have length equal to the maximum rank of the views provided
120  /*
121  Examples:
122  1. vector dot product in third dimension: {DIMENSION_MATCH, DIMENSION_MATCH, TENSOR_CONTRACTION}
123  - view1 and view2 should both be rank 3, and should match in all dimensions
124  - tensor_view should be rank 2, and should match view1 and view2 in first two dimensions
125  2. vector outer product in third dimension: {DIMENSION_MATCH, DIMENSION_MATCH, TENSOR_PRODUCT}
126  - view1 and view2 should both be rank 3, and should match in first two dimensions
127  - tensor_view should be rank 3, and should match view1 and view2 in first two dimensions
128  - in third dimension, tensor_view should have dimension equal to the product of the third dimension of view1 and the third dimension of view2
129  3. rank-3 view1 treated as vector times scalar rank-2 view2: {DIMENSION_MATCH, DIMENSION_MATCH, TENSOR_PRODUCT}
130  - here, the rank-2 view2 is interpreted as having an extent 1 third dimension
131 
132  We only allow TENSOR_CONTRACTION in final dimension(s)
133  */
134  // check that the above rules are satisfied:
135  unsigned max_component_rank = (view1.rank() > view2.rank()) ? view1.rank() : view2.rank();
136  unsigned max_rank = (tensor_view.rank() > max_component_rank) ? tensor_view.rank() : max_component_rank;
137 
138  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(rank_combination_types.extent(0) != max_rank, std::invalid_argument, "need to provide RankCombinationType for the largest-rank View");
139 
140  unsigned expected_rank = 0;
141  bool contracting = false;
142  for (unsigned d=0; d<rank_combination_types.extent(0); d++)
143  {
144  if (rank_combination_types[d] == TENSOR_CONTRACTION)
145  {
146  // check that view1 and view2 agree on the length of this dimension
147  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(view1.extent_int(d) != view2.extent_int(d), std::invalid_argument, "Contractions can only occur along ranks of equal length");
148  contracting = true;
149  }
150  else
151  {
152  if (!contracting)
153  {
154  expected_rank++;
155  if (rank_combination_types[d] == TENSOR_PRODUCT)
156  {
157  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(tensor_view.extent_int(d) != view1.extent_int(d) * view2.extent_int(d), std::invalid_argument, "For TENSOR_PRODUCT rank combination, the tensor View must have length in that dimension equal to the product of the two component views in that dimension");
158  }
159  else // matching
160  {
161  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(view1.extent_int(d) != view2.extent_int(d), std::invalid_argument, "For DIMENSION_MATCH rank combination, all three views must have length equal to each other in that rank");
162  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(tensor_view.extent_int(d) != view1.extent_int(d), std::invalid_argument, "For DIMENSION_MATCH rank combination, all three views must have length equal to each other in that rank");
163  }
164  }
165  else
166  {
167  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(contracting, std::invalid_argument, "encountered a non-contraction rank combination after a contraction; contractions can only go at the end");
168  }
169  }
170  }
171  INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(expected_rank != tensor_view.rank(), std::invalid_argument, "Tensor view does not match expected rank");
172  }
173 
176  KOKKOS_INLINE_FUNCTION
178  {
179  int view2_next_increment_rank = view2_iterator_.nextIncrementRank();
180  int view1_next_increment_rank = view1_iterator_.nextIncrementRank();
181  if (view2_next_increment_rank > view1_next_increment_rank) return view2_next_increment_rank;
182  else return view1_next_increment_rank;
183  }
184 
187  KOKKOS_INLINE_FUNCTION
188  int increment()
189  {
190  // proceed to the next view1/view2 combination
191  // where we're doing a dimension match, then all three iterators should increment in tandem
192  // where we're doing a contraction, view1/view2 should increment in tandem, while tensor_view should be fixed
193  // where we're doing a tensor product, view1 and tensor_view increment in tandem, while view2 is fixed
194 
195  // note that regardless of the choice, view1 should be incremented, with one exception:
196  // If we are doing a tensor product, then view1 can be understood to be in an interior for loop, and it should loop around.
197  // We can detect this by checking which the least rank that would be updated -- if view2's least rank exceeds view1's, then:
198  // - view1 should be reset, AND
199  // - view2 should be incremented (as should the tensor view)
200  int view2_next_increment_rank = view2_iterator_.nextIncrementRank();
201  int view1_next_increment_rank = view1_iterator_.nextIncrementRank();
202  if (view2_next_increment_rank > view1_next_increment_rank)
203  {
204  // if we get here, we should be doing a tensor product in the view2 rank that will change
205  device_assert(rank_combination_types_[view2_next_increment_rank]==TENSOR_PRODUCT);
206  view1_iterator_.reset(view2_next_increment_rank); // set to 0 from the tensor product rank inward -- this is "looping around"
207  view2_iterator_.increment();
208  tensor_view_iterator_.increment();
209  return view2_next_increment_rank;
210  }
211  else
212  {
213  int view1_rank_change = view1_iterator_.increment();
214  if (view1_rank_change >= 0)
215  {
216  switch (rank_combination_types_[view1_rank_change])
217  {
218  case DIMENSION_MATCH:
219  view2_iterator_.increment();
220  tensor_view_iterator_.increment();
221  break;
222  case TENSOR_PRODUCT:
223  // view1 increments fastest; the only time we increment view2 is when view1 loops around; we handle that above
224  tensor_view_iterator_.increment();
225  break;
226  case TENSOR_CONTRACTION:
227  // view1 and view2 increment in tandem; we don't increment tensor_view while contraction is taking place
228  view2_iterator_.increment();
229  }
230  }
231  return view1_rank_change;
232  }
233  }
234 
237  KOKKOS_INLINE_FUNCTION
238  void setLocation(const Kokkos::Array<int,7> location)
239  {
240  view1_iterator_.setLocation(location);
241  view2_iterator_.setLocation(location);
242  tensor_view_iterator_.setLocation(location);
243  }
244 
248  KOKKOS_INLINE_FUNCTION
249  void setLocation(Kokkos::Array<int,7> location1, Kokkos::Array<int,7> location2)
250  {
251  view1_iterator_.setLocation(location1);
252  view2_iterator_.setLocation(location2);
253  Kokkos::Array<int,7> tensor_location = location1;
254  for (unsigned d=0; d<rank_combination_types_.extent(0); d++)
255  {
256  switch (rank_combination_types_[d])
257  {
258  case TENSOR_PRODUCT:
259  // view1 index is fastest-moving:
260  tensor_location[d] = location2[d] * view1_iterator_.getExtent(d) + location1[d];
261  break;
262  case DIMENSION_MATCH:
263  // we copied location1 into tensor_location to initialize -- that's correct in this dimension
264  break;
265  case TENSOR_CONTRACTION:
266  tensor_location[d] = 0;
267  break;
268  }
269  }
270 #ifdef HAVE_INTREPID2_DEBUG
271  // check that the location makes sense
272  for (unsigned d=0; d<rank_combination_types_.extent(0); d++)
273  {
274  switch (rank_combination_types_[d])
275  {
276  case TENSOR_PRODUCT:
277  // in this case, the two indices are independent
278  break;
279  case DIMENSION_MATCH:
280  case TENSOR_CONTRACTION:
281  device_assert(location1[d] == location2[d]);
282  break;
283  }
284  // let's check that the indices are in bounds:
285  device_assert(location1[d] < view1_iterator_.getExtent(d));
286  device_assert(location2[d] < view2_iterator_.getExtent(d));
287  device_assert(tensor_location[d] < tensor_view_iterator_.getExtent(d));
288  }
289 #endif
290  tensor_view_iterator_.setLocation(tensor_location);
291  }
292 
295  KOKKOS_INLINE_FUNCTION
296  ScalarType getView1Entry()
297  {
298  return view1_iterator_.get();
299  }
300 
303  KOKKOS_INLINE_FUNCTION
304  ScalarType getView2Entry()
305  {
306  return view2_iterator_.get();
307  }
308 
311  KOKKOS_INLINE_FUNCTION
312  void set(ScalarType value)
313  {
314  tensor_view_iterator_.set(value);
315  }
316  };
317 
318 } // namespace Intrepid2
319 
320 #endif /* Intrepid2_TensorViewIterator_h */
KOKKOS_INLINE_FUNCTION ScalarType getView1Entry()
KOKKOS_INLINE_FUNCTION void reset(unsigned from_rank_number=0)
KOKKOS_INLINE_FUNCTION int nextIncrementRank()
KOKKOS_INLINE_FUNCTION int getExtent(int dimension)
KOKKOS_INLINE_FUNCTION TensorViewIterator(TensorViewType tensor_view, ViewType1 view1, ViewType2 view2, RankCombinationViewType rank_combination_types)
Constructor.
KOKKOS_INLINE_FUNCTION ScalarType get()
KOKKOS_INLINE_FUNCTION void setLocation(Kokkos::Array< int, 7 > location1, Kokkos::Array< int, 7 > location2)
Implementation of an assert that can safely be called from device code.
Iterator allows linear traversal of (part of) a Kokkos View in a manner that is agnostic to its rank...
KOKKOS_INLINE_FUNCTION void set(const ScalarType &value)
KOKKOS_INLINE_FUNCTION void setLocation(const Kokkos::Array< int, 7 > location)
KOKKOS_INLINE_FUNCTION ScalarType getView2Entry()
KOKKOS_INLINE_FUNCTION int increment()
KOKKOS_INLINE_FUNCTION void setLocation(const Kokkos::Array< int, 7 > &location)
KOKKOS_INLINE_FUNCTION int nextIncrementRank()
A helper class that allows iteration over three Kokkos Views simultaneously, according to tensor comb...