Intrepid2
Intrepid2_ArrayToolsDefScalar.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Intrepid2 Package
4//
5// Copyright 2007 NTESS and the Intrepid2 contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
16#ifndef __INTREPID2_ARRAYTOOLS_DEF_SCALAR_HPP__
17#define __INTREPID2_ARRAYTOOLS_DEF_SCALAR_HPP__
18
19namespace Intrepid2 {
20
21 namespace FunctorArrayTools {
25 template <typename OutputViewType,
26 typename inputLeftViewType,
27 typename inputRightViewType,
28 bool equalRank,
29 bool reciprocal>
31 {
32 OutputViewType _output;
33 const inputLeftViewType _inputLeft;
34 const inputRightViewType _inputRight;
35
36 KOKKOS_INLINE_FUNCTION
37 F_scalarMultiplyScalar(OutputViewType output_,
38 inputLeftViewType inputLeft_,
39 inputRightViewType inputRight_)
40 : _output(output_),
41 _inputLeft(inputLeft_),
42 _inputRight(inputRight_) {}
43
44 // DataData
45 KOKKOS_INLINE_FUNCTION
46 void operator()(const ordinal_type cl, const ordinal_type pt) const
47 {
48 const auto val = _inputLeft(cl, pt % _inputLeft.extent(1));
49
50 if constexpr (equalRank) {
51 if constexpr (reciprocal)
52 _output(cl, pt) = _inputRight(cl, pt) / val;
53 else
54 _output(cl, pt) = val * _inputRight(cl, pt);
55 }
56 else {
57 if constexpr (reciprocal)
58 _output(cl, pt) = _inputRight(pt) / val;
59 else
60 _output(cl, pt) = val * _inputRight(pt);
61 }
62 }
63
64 // DataField
65 KOKKOS_INLINE_FUNCTION
66 void operator()(const ordinal_type cl, const ordinal_type bf, const ordinal_type pt) const
67 {
68 const auto val = _inputLeft(cl, pt % _inputLeft.extent(1));
69
70 if constexpr (equalRank) {
71 if constexpr (reciprocal)
72 _output(cl, bf, pt) = _inputRight(cl, bf, pt) / val;
73 else
74 _output(cl, bf, pt) = val * _inputRight(cl, bf, pt);
75 }
76 else {
77 if constexpr (reciprocal)
78 _output(cl, bf, pt) = _inputRight(bf, pt) / val;
79 else
80 _output(cl, bf, pt) = val * _inputRight(bf, pt);
81 }
82 }
83 };
84
85 template <typename OutputViewType,
86 typename inputLeftViewType,
87 typename inputRightViewType,
88 bool equalRank,
89 bool reciprocal>
91 {
92 OutputViewType _output;
93 const inputLeftViewType _inputLeft;
94 const inputRightViewType _inputRight;
95
96 KOKKOS_INLINE_FUNCTION
97 F_scalarMultiplyVector(OutputViewType output_,
98 inputLeftViewType inputLeft_,
99 inputRightViewType inputRight_)
100 : _output(output_),
101 _inputLeft(inputLeft_),
102 _inputRight(inputRight_) {}
103
104 // DataData
105 KOKKOS_INLINE_FUNCTION
106 void operator()(const ordinal_type cl, const ordinal_type pt) const
107 {
108 const auto val = _inputLeft(cl, pt % _inputLeft.extent(1));
109
110 if constexpr (equalRank) {
111 if constexpr (reciprocal) {
112 for (ordinal_type i = 0; i < _output.extent_int(2); ++i)
113 _output(cl, pt, i) = _inputRight(cl, pt, i) / val;
114 }
115 else {
116 for (ordinal_type i = 0; i < _output.extent_int(2); ++i)
117 _output(cl, pt, i) = val * _inputRight(cl, pt, i);
118 }
119 }
120 else {
121 if constexpr (reciprocal) {
122 for (ordinal_type i = 0; i < _output.extent_int(2); ++i)
123 _output(cl, pt, i) = _inputRight(pt, i) / val;
124 }
125 else {
126 for (ordinal_type i = 0; i < _output.extent_int(2); ++i)
127 _output(cl, pt, i) = val * _inputRight(pt, i);
128 }
129 }
130 }
131
132 // DataField
133 KOKKOS_INLINE_FUNCTION
134 void operator()(const ordinal_type cl, const ordinal_type bf, const ordinal_type pt) const
135 {
136 const auto val = _inputLeft(cl, pt % _inputLeft.extent(1));
137
138 if constexpr (equalRank) {
139 if constexpr (reciprocal) {
140 for (ordinal_type i = 0; i < _output.extent_int(3); ++i)
141 _output(cl, bf, pt, i) = _inputRight(cl, bf, pt, i) / val;
142 }
143 else {
144 for (ordinal_type i = 0; i < _output.extent_int(3); ++i)
145 _output(cl, bf, pt, i) = val * _inputRight(cl, bf, pt, i);
146 }
147 }
148 else {
149 if constexpr (reciprocal) {
150 for (ordinal_type i = 0; i < _output.extent_int(3); ++i)
151 _output(cl, bf, pt, i) = _inputRight(bf, pt, i) / val;
152 }
153 else {
154 for (ordinal_type i = 0; i < _output.extent_int(3); ++i)
155 _output(cl, bf, pt, i) = val * _inputRight(bf, pt, i);
156 }
157 }
158 }
159 };
160
161 template <typename OutputViewType,
162 typename inputLeftViewType,
163 typename inputRightViewType,
164 bool equalRank,
165 bool reciprocal>
167 {
168 OutputViewType _output;
169 const inputLeftViewType _inputLeft;
170 const inputRightViewType _inputRight;
171
172 KOKKOS_INLINE_FUNCTION
173 F_scalarMultiplyTensor(OutputViewType output_,
174 inputLeftViewType inputLeft_,
175 inputRightViewType inputRight_)
176 : _output(output_),
177 _inputLeft(inputLeft_),
178 _inputRight(inputRight_) {}
179
180 // DataData
181 KOKKOS_INLINE_FUNCTION
182 void operator()(const ordinal_type cl, const ordinal_type pt) const
183 {
184 const auto val = _inputLeft(cl, pt % _inputLeft.extent(1));
185
186 if constexpr (equalRank) {
187 if constexpr (reciprocal) {
188 for (ordinal_type i = 0; i < _output.extent_int(2); ++i)
189 for (ordinal_type j = 0; j < _output.extent_int(3); ++j)
190 _output(cl, pt, i, j) = _inputRight(cl, pt, i, j) / val;
191 }
192 else {
193 for (ordinal_type i = 0; i < _output.extent_int(2); ++i)
194 for (ordinal_type j = 0; j < _output.extent_int(3); ++j)
195 _output(cl, pt, i, j) = val * _inputRight(cl, pt, i, j);
196 }
197 }
198 else {
199 if constexpr (reciprocal) {
200 for (ordinal_type i = 0; i < _output.extent_int(2); ++i)
201 for (ordinal_type j = 0; j < _output.extent_int(3); ++j)
202 _output(cl, pt, i, j) = _inputRight(pt, i, j) / val;
203 }
204 else {
205 for (ordinal_type i = 0; i < _output.extent_int(2); ++i)
206 for (ordinal_type j = 0; j < _output.extent_int(3); ++j)
207 _output(cl, pt, i, j) = val * _inputRight(pt, i, j);
208 }
209 }
210 }
211
212 // DataField
213 KOKKOS_INLINE_FUNCTION
214 void operator()(const ordinal_type cl, const ordinal_type bf, const ordinal_type pt) const
215 {
216 const auto val = _inputLeft(cl, pt % _inputLeft.extent(1));
217
218 if constexpr (equalRank) {
219 if constexpr (reciprocal) {
220 for (ordinal_type i = 0; i < _output.extent_int(3); ++i)
221 for (ordinal_type j = 0; j < _output.extent_int(4); ++j)
222 _output(cl, bf, pt, i, j) = _inputRight(cl, bf, pt, i, j) / val;
223 }
224 else {
225 for (ordinal_type i = 0; i < _output.extent_int(3); ++i)
226 for (ordinal_type j = 0; j < _output.extent_int(4); ++j)
227 _output(cl, bf, pt, i, j) = val * _inputRight(cl, bf, pt, i, j);
228 }
229 }
230 else {
231 if constexpr (reciprocal) {
232 for (ordinal_type i = 0; i < _output.extent_int(3); ++i)
233 for (ordinal_type j = 0; j < _output.extent_int(4); ++j)
234 _output(cl, bf, pt, i, j) = _inputRight(bf, pt, i, j) / val;
235 }
236 else {
237 for (ordinal_type i = 0; i < _output.extent_int(3); ++i)
238 for (ordinal_type j = 0; j < _output.extent_int(4); ++j)
239 _output(cl, bf, pt, i, j) = val * _inputRight(bf, pt, i, j);
240 }
241 }
242 }
243 };
244 }
245
246 template<typename DeviceType>
247 template<typename outputFieldValueType, class ...outputFieldProperties,
248 typename inputDataValueType, class ...inputDataProperties,
249 typename inputFieldValueType, class ...inputFieldProperties>
250 void
252 scalarMultiplyDataField( Kokkos::DynRankView<outputFieldValueType,outputFieldProperties...> outputFields,
253 const Kokkos::DynRankView<inputDataValueType, inputDataProperties...> inputData,
254 const Kokkos::DynRankView<inputFieldValueType, inputFieldProperties...> inputFields,
255 const bool reciprocal ) {
256
257#ifdef HAVE_INTREPID2_DEBUG
258 {
259 INTREPID2_TEST_FOR_EXCEPTION( inputData.rank() != 2, std::invalid_argument,
260 ">>> ERROR (ArrayTools::scalarMultiplyDataField): Input data container must have rank 2.");
261
262 if (outputFields.rank() <= inputFields.rank()) {
263 INTREPID2_TEST_FOR_EXCEPTION( inputFields.rank() < 3 || inputFields.rank() > 5, std::invalid_argument,
264 ">>> ERROR (ArrayTools::scalarMultiplyDataField): Input fields container must have rank 3, 4, or 5.");
265 INTREPID2_TEST_FOR_EXCEPTION( outputFields.rank() != inputFields.rank(), std::invalid_argument,
266 ">>> ERROR (ArrayTools::scalarMultiplyDataField): Input and output fields containers must have the same rank.");
267 INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(0) != inputData.extent(0), std::invalid_argument,
268 ">>> ERROR (ArrayTools::scalarMultiplyDataField): Zeroth dimensions (number of integration domains) of the fields and data input containers must agree!");
269 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(1) != inputFields.extent(2) &&
270 inputData.extent(1) != 1, std::invalid_argument,
271 ">>> ERROR (ArrayTools::scalarMultiplyDataField): Second dimension of the fields input container and first dimension of data input container (number of integration points) must agree, or first data dimension must be 1!");
272 for (size_type i=0;i<inputFields.rank();++i) {
273 INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(i) != outputFields.extent(i), std::invalid_argument,
274 ">>> ERROR (ArrayTools::scalarMultiplyDataField): inputFields dimension (i) does not match to the dimension (i) of outputFields");
275 }
276 }
277 else {
278 INTREPID2_TEST_FOR_EXCEPTION( inputFields.rank() < 2 || inputFields.rank() > 4, std::invalid_argument,
279 ">>> ERROR (ArrayTools::scalarMultiplyDataField): Input fields container must have rank 2, 3, or 4.");
280 INTREPID2_TEST_FOR_EXCEPTION( outputFields.rank() != (inputFields.rank()+1), std::invalid_argument,
281 ">>> ERROR (ArrayTools::scalarMultiplyDataField): The rank of the input fields container must be one less than the rank of the output fields container.");
282 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(1) != inputFields.extent(1) &&
283 inputData.extent(1) != 1, std::invalid_argument,
284 ">>> ERROR (ArrayTools::scalarMultiplyDataField): First dimensions of fields input container and data input container (number of integration points) must agree or first data dimension must be 1!");
285 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(0) != outputFields.extent(0), std::invalid_argument,
286 ">>> ERROR (ArrayTools::scalarMultiplyDataField): Zeroth dimensions of fields output container and data input containers (number of integration domains) must agree!");
287 for (size_type i=0;i<inputFields.rank();++i) {
288 INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(i) != outputFields.extent(i+1), std::invalid_argument,
289 ">>> ERROR (ArrayTools::scalarMultiplyDataField): inputFields dimension (i) does not match to the dimension (i+1) of outputFields");
290 }
291 }
292 }
293#endif
294
295 typedef Kokkos::DynRankView<outputFieldValueType,outputFieldProperties...> outputFieldViewType;
296 typedef Kokkos::DynRankView<inputDataValueType,inputDataProperties...> inputDataViewType;
297 typedef Kokkos::DynRankView<inputFieldValueType,inputFieldProperties...> inputFieldViewType;
298
299 using range_policy_type = Kokkos::MDRangePolicy
300 < ExecSpaceType, Kokkos::Rank<3>, Kokkos::IndexType<ordinal_type> >;
301
302 const range_policy_type policy( { 0, 0, 0 },
303 { /*C*/ outputFields.extent(0), /*F*/ outputFields.extent(1), /*P*/ outputFields.extent(2) } );
304
305 const bool equalRank = ( outputFields.rank() == inputFields.rank() );
306
307 if (outputFields.rank() == 3) {
308 if (equalRank) {
309 if (reciprocal) {
311 Kokkos::parallel_for( policy, FunctorType(outputFields, inputData, inputFields) );
312 } else {
314 Kokkos::parallel_for( policy, FunctorType(outputFields, inputData, inputFields) );
315 }
316 }
317 else {
318 if (reciprocal) {
320 Kokkos::parallel_for( policy, FunctorType(outputFields, inputData, inputFields) );
321 } else {
323 Kokkos::parallel_for( policy, FunctorType(outputFields, inputData, inputFields) );
324 }
325 }
326 }
327 else if (outputFields.rank() == 4) {
328 if (equalRank) {
329 if (reciprocal) {
331 Kokkos::parallel_for( policy, FunctorType(outputFields, inputData, inputFields) );
332 } else {
334 Kokkos::parallel_for( policy, FunctorType(outputFields, inputData, inputFields) );
335 }
336 }
337 else {
338 if (reciprocal) {
340 Kokkos::parallel_for( policy, FunctorType(outputFields, inputData, inputFields) );
341 } else {
343 Kokkos::parallel_for( policy, FunctorType(outputFields, inputData, inputFields) );
344 }
345 }
346 }
347 else {
348 if (equalRank) {
349 if (reciprocal) {
351 Kokkos::parallel_for( policy, FunctorType(outputFields, inputData, inputFields) );
352 } else {
354 Kokkos::parallel_for( policy, FunctorType(outputFields, inputData, inputFields) );
355 }
356 }
357 else {
358 if (reciprocal) {
360 Kokkos::parallel_for( policy, FunctorType(outputFields, inputData, inputFields) );
361 } else {
363 Kokkos::parallel_for( policy, FunctorType(outputFields, inputData, inputFields) );
364 }
365 }
366 }
367 }
368
369
370 template<typename DeviceType>
371 template<typename outputDataValueType, class ...outputDataProperties,
372 typename inputDataLeftValueType, class ...inputDataLeftProperties,
373 typename inputDataRightValueType, class ...inputDataRightProperties>
374 void
376 scalarMultiplyDataData( Kokkos::DynRankView<outputDataValueType, outputDataProperties...> outputData,
377 const Kokkos::DynRankView<inputDataLeftValueType, inputDataLeftProperties...> inputDataLeft,
378 const Kokkos::DynRankView<inputDataRightValueType,inputDataRightProperties...> inputDataRight,
379 const bool reciprocal ) {
380
381#ifdef HAVE_INTREPID2_DEBUG
382 {
383 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.rank() != 2, std::invalid_argument,
384 ">>> ERROR (ArrayTools::scalarMultiplyDataData): Left input data container must have rank 2.");
385
386 if (outputData.rank() <= inputDataRight.rank()) {
387 INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.rank() < 2 || inputDataRight.rank() > 4, std::invalid_argument,
388 ">>> ERROR (ArrayTools::scalarMultiplyDataData): Right input data container must have rank 2, 3, or 4.");
389 INTREPID2_TEST_FOR_EXCEPTION( outputData.rank() != inputDataRight.rank(), std::invalid_argument,
390 ">>> ERROR (ArrayTools::scalarMultiplyDataData): Right input and output data containers must have the same rank.");
391 INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.extent(0) != inputDataLeft.extent(0), std::invalid_argument,
392 ">>> ERROR (ArrayTools::scalarMultiplyDataData): Zeroth dimensions (number of integration domains) of the left and right data input containers must agree!");
393 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(1) != inputDataRight.extent(1) &&
394 inputDataLeft.extent(1) != 1, std::invalid_argument,
395 ">>> ERROR (ArrayTools::scalarMultiplyDataData): First dimensions of the left and right data input containers (number of integration points) must agree, or first dimension of the left data input container must be 1!");
396 for (size_type i=0;i<inputDataRight.rank();++i) {
397 INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.extent(i) != outputData.extent(i), std::invalid_argument,
398 ">>> ERROR (ArrayTools::scalarMultiplyDataData): inputDataRight dimension (i) does not match to the dimension (i) of outputData");
399 }
400 } else {
401 INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.rank() < 1 || inputDataRight.rank() > 3, std::invalid_argument,
402 ">>> ERROR (ArrayTools::scalarMultiplyDataData): Right input data container must have rank 1, 2, or 3.");
403 INTREPID2_TEST_FOR_EXCEPTION( outputData.rank() != (inputDataRight.rank()+1), std::invalid_argument,
404 ">>> ERROR (ArrayTools::scalarMultiplyDataData): The rank of the right input data container must be one less than the rank of the output data container.");
405 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(1) != inputDataRight.extent(0) &&
406 inputDataLeft.extent(1) != 1, std::invalid_argument,
407 ">>> ERROR (ArrayTools::scalarMultiplyDataData): Zeroth dimension of the right input data container and first dimension of the left data input container (number of integration points) must agree or first dimension of the left data input container must be 1!");
408 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(0) != outputData.extent(0), std::invalid_argument,
409 ">>> ERROR (ArrayTools::scalarMultiplyDataData): Zeroth dimensions of data output and left data input containers (number of integration domains) must agree!");
410 for (size_type i=0;i<inputDataRight.rank();++i) {
411 INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.extent(i) != outputData.extent(i+1), std::invalid_argument,
412 ">>> ERROR (ArrayTools::scalarMultiplyDataData): inputDataRight dimension (i) does not match to the dimension (i+1) of outputData");
413 }
414 }
415 }
416#endif
417
418 typedef Kokkos::DynRankView<outputDataValueType,outputDataProperties...> outputDataViewType;
419 typedef Kokkos::DynRankView<inputDataLeftValueType,inputDataLeftProperties...> inputDataLeftViewType;
420 typedef Kokkos::DynRankView<inputDataRightValueType,inputDataRightProperties...> inputDataRightViewType;
421
422 using range_policy_type = Kokkos::MDRangePolicy
423 < ExecSpaceType, Kokkos::Rank<2>, Kokkos::IndexType<ordinal_type> >;
424
425 const range_policy_type policy( { 0, 0 },
426 { /*C*/ outputData.extent(0), /*P*/ outputData.extent(1) } );
427
428 const bool equalRank = ( outputData.rank() == inputDataRight.rank() );
429 if (outputData.rank() == 2) {
430 if (equalRank) {
431 if (reciprocal) {
433 Kokkos::parallel_for( policy, FunctorType(outputData, inputDataLeft, inputDataRight) );
434 } else {
436 Kokkos::parallel_for( policy, FunctorType(outputData, inputDataLeft, inputDataRight) );
437 }
438 }
439 else {
440 if (reciprocal) {
442 Kokkos::parallel_for( policy, FunctorType(outputData, inputDataLeft, inputDataRight) );
443 } else {
445 Kokkos::parallel_for( policy, FunctorType(outputData, inputDataLeft, inputDataRight) );
446 }
447 }
448 }
449 else if (outputData.rank() == 3) {
450 if (equalRank) {
451 if (reciprocal) {
453 Kokkos::parallel_for( policy, FunctorType(outputData, inputDataLeft, inputDataRight) );
454 } else {
456 Kokkos::parallel_for( policy, FunctorType(outputData, inputDataLeft, inputDataRight) );
457 }
458 }
459 else {
460 if (reciprocal) {
462 Kokkos::parallel_for( policy, FunctorType(outputData, inputDataLeft, inputDataRight) );
463 } else {
465 Kokkos::parallel_for( policy, FunctorType(outputData, inputDataLeft, inputDataRight) );
466 }
467 }
468 }
469 else {
470 if (equalRank) {
471 if (reciprocal) {
473 Kokkos::parallel_for( policy, FunctorType(outputData, inputDataLeft, inputDataRight) );
474 } else {
476 Kokkos::parallel_for( policy, FunctorType(outputData, inputDataLeft, inputDataRight) );
477 }
478 }
479 else {
480 if (reciprocal) {
482 Kokkos::parallel_for( policy, FunctorType(outputData, inputDataLeft, inputDataRight) );
483 } else {
485 Kokkos::parallel_for( policy, FunctorType(outputData, inputDataLeft, inputDataRight) );
486 }
487 }
488 }
489 }
490
491} // end namespace Intrepid2
492
493#endif
static void scalarMultiplyDataField(Kokkos::DynRankView< outputFieldValueType, outputFieldProperties... > outputFields, const Kokkos::DynRankView< inputDataValueType, inputDataProperties... > inputData, const Kokkos::DynRankView< inputFieldValueType, inputFieldProperties... > inputFields, const bool reciprocal=false)
There are two use cases: (1) multiplies a rank-3, 4, or 5 container inputFields with dimensions (C,...
static void scalarMultiplyDataData(Kokkos::DynRankView< outputDataValueType, outputDataProperties... > outputData, const Kokkos::DynRankView< inputDataLeftValueType, inputDataLeftProperties... > inputDataLeft, const Kokkos::DynRankView< inputDataRightValueType, inputDataRightProperties... > inputDataRight, const bool reciprocal=false)
There are two use cases: (1) multiplies a rank-2, 3, or 4 container inputDataRight with dimensions (C...
Functors for scalarMultiply see Intrepid2::ArrayTools for more.