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
58namespace 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
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 */
Implementation of an assert that can safely be called from device code.
KOKKOS_INLINE_FUNCTION void device_assert(bool val)
#define INTREPID2_TEST_FOR_EXCEPTION_DEVICE_SAFE(test, x, msg)
Iterator allows linear traversal of (part of) a Kokkos View in a manner that is agnostic to its rank.
A helper class that allows iteration over three Kokkos Views simultaneously, according to tensor comb...
KOKKOS_INLINE_FUNCTION void setLocation(Kokkos::Array< int, 7 > location1, Kokkos::Array< int, 7 > location2)
KOKKOS_INLINE_FUNCTION ScalarType getView1Entry()
KOKKOS_INLINE_FUNCTION int nextIncrementRank()
KOKKOS_INLINE_FUNCTION TensorViewIterator(TensorViewType tensor_view, ViewType1 view1, ViewType2 view2, RankCombinationViewType rank_combination_types)
Constructor.
KOKKOS_INLINE_FUNCTION void setLocation(const Kokkos::Array< int, 7 > location)
KOKKOS_INLINE_FUNCTION void set(ScalarType value)
KOKKOS_INLINE_FUNCTION ScalarType getView2Entry()
A helper class that allows iteration over some part of a Kokkos View, while allowing the calling code...
KOKKOS_INLINE_FUNCTION int nextIncrementRank()
KOKKOS_INLINE_FUNCTION int getExtent(int dimension)
KOKKOS_INLINE_FUNCTION void set(const ScalarType &value)
KOKKOS_INLINE_FUNCTION int increment()
KOKKOS_INLINE_FUNCTION ScalarType get()
KOKKOS_INLINE_FUNCTION void setLocation(const Kokkos::Array< int, 7 > &location)
KOKKOS_INLINE_FUNCTION void reset(unsigned from_rank_number=0)