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 
19 namespace 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
Functors for scalarMultiply see Intrepid2::ArrayTools for more.
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...
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...