Intrepid2
Intrepid2_ArrayToolsDefTensor.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_TENSOR_HPP__
17 #define __INTREPID2_ARRAYTOOLS_DEF_TENSOR_HPP__
18 
19 namespace Intrepid2 {
20 
21  namespace FunctorArrayTools {
25  template < typename OutputViewType, typename leftInputViewType, typename rightInputViewType >
26  struct F_crossProduct{
27  OutputViewType _output;
28  const leftInputViewType _leftInput;
29  const rightInputViewType _rightInput;
30  const bool _hasField, _isCrossProd3D;
31 
32  KOKKOS_INLINE_FUNCTION
33  F_crossProduct(OutputViewType output_,
34  leftInputViewType leftInput_,
35  rightInputViewType rightInput_,
36  const bool hasField_,
37  const bool isCrossProd3D_)
38  : _output(output_), _leftInput(leftInput_), _rightInput(rightInput_),
39  _hasField(hasField_), _isCrossProd3D(isCrossProd3D_) {}
40 
41  KOKKOS_INLINE_FUNCTION
42  void operator()(const size_type iter) const {
43  size_type cell, field(0), point;
44  size_type rightRank = _rightInput.rank();
45 
46  if (_hasField)
47  unrollIndex( cell, field, point,
48  _output.extent(0),
49  _output.extent(1),
50  _output.extent(2),
51  iter );
52  else
53  unrollIndex( cell, point,
54  _output.extent(0),
55  _output.extent(1),
56  iter );
57 
58  auto result = ( _hasField ? Kokkos::subview(_output, cell, field, point, Kokkos::ALL()) :
59  Kokkos::subview(_output, cell, point, Kokkos::ALL()));
60 
61  auto left = Kokkos::subview(_leftInput, cell, point, Kokkos::ALL());
62 
63  auto right = (rightRank == 2 + size_type(_hasField)) ?
64  ( _hasField ? Kokkos::subview(_rightInput, field, point, Kokkos::ALL()) :
65  Kokkos::subview(_rightInput, point, Kokkos::ALL())) :
66  ( _hasField ? Kokkos::subview(_rightInput, cell, field, point, Kokkos::ALL()) :
67  Kokkos::subview(_rightInput, cell, point, Kokkos::ALL()));
68 
69  // branch prediction is possible
70  if (_isCrossProd3D) {
71  result(0) = left(1)*right(2) - left(2)*right(1);
72  result(1) = left(2)*right(0) - left(0)*right(2);
73  result(2) = left(0)*right(1) - left(1)*right(0);
74  } else {
75  result(0) = left(0)*right(1) - left(1)*right(0);
76  }
77  }
78  };
79  } //namespace
80 
81  template<typename DeviceType>
82  template<typename outputValueType, class ...outputProperties,
83  typename leftInputValueType, class ...leftInputProperties,
84  typename rightInputValueType, class ...rightInputProperties>
85  void
86  ArrayTools<DeviceType>::Internal::
87  crossProduct( Kokkos::DynRankView<outputValueType, outputProperties...> output,
88  const Kokkos::DynRankView<leftInputValueType, leftInputProperties...> leftInput,
89  const Kokkos::DynRankView<rightInputValueType,rightInputProperties...> rightInput,
90  const bool hasField ) {
91 
92  typedef Kokkos::DynRankView<outputValueType, outputProperties...> OutputViewType;
93  typedef const Kokkos::DynRankView<leftInputValueType, leftInputProperties...> leftInputViewType;
94  typedef const Kokkos::DynRankView<rightInputValueType,rightInputProperties...> rightInputViewType;
95  typedef FunctorArrayTools::F_crossProduct<OutputViewType, leftInputViewType, rightInputViewType> FunctorType;
96 
97  const size_type loopSize = ( hasField ? output.extent(0)*output.extent(1)*output.extent(2) :
98  output.extent(0)*output.extent(1) );
99  const bool isCrossProd3D = (leftInput.extent(2) == 3);
100  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
101  Kokkos::parallel_for( policy, FunctorType(output, leftInput, rightInput, hasField, isCrossProd3D) );
102  }
103 
104 
105 
106  template<typename DeviceType>
107  template<typename outputFieldValueType, class ...outputFieldProperties,
108  typename inputDataValueType, class ...inputDataProperties,
109  typename inputFieldValueType, class ...inputFieldProperties>
110  void
111  ArrayTools<DeviceType>::
112  crossProductDataField( Kokkos::DynRankView<outputFieldValueType,outputFieldProperties...> outputFields,
113  const Kokkos::DynRankView<inputDataValueType, inputDataProperties...> inputData,
114  const Kokkos::DynRankView<inputFieldValueType, inputFieldProperties...> inputFields ) {
115 
116 #ifdef HAVE_INTREPID2_DEBUG
117  {
118  /*
119  * Check array rank and spatial dimension range (if applicable)
120  * (1) inputData(C,P,D);
121  * (2) inputFields(C,F,P,D) or (F,P,D);
122  * (3) outputFields(C,F,P,D) in 3D, or (C,F,P) in 2D
123  */
124  // (1) inputData is (C, P, D) and 2 <= D <= 3 is required
125  INTREPID2_TEST_FOR_EXCEPTION( inputData.rank() != 3, std::invalid_argument,
126  ">>> ERROR (ArrayTools::crossProductDataField): inputData must have rank 3");
127  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(2) < 2 || inputData.extent(2) > 3, std::invalid_argument,
128  ">>> ERROR (ArrayTools::crossProductDataField): inputData dimension(2) must be 2 or 3");
129 
130  // (2) inputFields is (C, F, P, D) or (F, P, D) and 2 <= (D=dimension(rank - 1)) <= 3 is required.
131  INTREPID2_TEST_FOR_EXCEPTION( inputFields.rank() < 3 || inputFields.rank() > 4, std::invalid_argument,
132  ">>> ERROR (ArrayTools::crossProductDataField): inputFields must have rank 3 or 4" );
133  INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(inputFields.rank()-1) < 2 ||
134  inputFields.extent(inputFields.rank()-1) > 3, std::invalid_argument,
135  ">>> ERROR (ArrayTools::crossProductDataField): inputFields dimension (rank-1) must have rank 2 or 3" );
136 
137  // (3) outputFields is (C,F,P,D) in 3D and (C,F,P) in 2D => rank = inputData.extent(2) + 1
138  INTREPID2_TEST_FOR_EXCEPTION( outputFields.rank() != inputData.extent(2)+1, std::invalid_argument,
139  ">>> ERROR (ArrayTools::crossProductDataField): outputFields rank must match to inputData dimension(2)+1");
140  /*
141  * Dimension cross-checks:
142  * (1) inputData vs. inputFields
143  * (2) outputFields vs. inputData
144  * (3) outputFields vs. inputFields
145  *
146  * Cross-check (1):
147  */
148  if (inputFields.rank() == 4) {
149  // inputData(C,P,D) vs. inputFields(C,F,P,D): dimensions C, P, D must match
150  const size_type f1[] = { 0, 1, 2 }, f2[] = { 0, 2, 3 };
151  for (size_type i=0; i<3; ++i) {
152  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
153  ">>> ERROR (ArrayTools::crossProductDataField): inputData dimension does not match with inputFields");
154  }
155  } else {
156  // inputData(C,P,D) vs. inputFields(F,P,D): dimensions P, D must match
157  const size_type f1[] = { 1, 2 }, f2[] = { 1, 2 };
158  for (size_type i=0; i<2; ++i) {
159  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
160  ">>> ERROR (ArrayTools::crossProductDataField): inputData dimension does not match with inputFields");
161  }
162  }
163  /*
164  * Cross-check (2):
165  */
166  if (inputData.extent(2) == 2) {
167  // in 2D: outputFields(C,F,P) vs. inputData(C,P,D): dimensions C,P must match
168  // inputData(C,P,D) vs. inputFields(F,P,D): dimensions P, D must match
169  const size_type f1[] = { 0, 2 }, f2[] = { 0, 1 };
170  for (size_type i=0; i<2; ++i) {
171  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputData.extent(f2[i]), std::invalid_argument,
172  ">>> ERROR (ArrayTools::crossProductDataField): outputFields dimension does not match with inputData");
173  }
174  } else {
175  // in 3D: outputFields(C,F,P,D) vs. inputData(C,P,D): dimensions C,P,D must match
176  const size_type f1[] = { 0, 2, 3 }, f2[] = { 0, 1, 2 };
177  for (size_type i=0; i<3; ++i) {
178  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputData.extent(f2[i]), std::invalid_argument,
179  ">>> ERROR (ArrayTools::crossProductDataField): outputFields dimension does not match with inputData");
180  }
181  }
182  /*
183  * Cross-check (3):
184  */
185  if (inputData.extent(2) == 2) {
186  // In 2D:
187  if (inputFields.rank() == 4) {
188  // and rank-4 inputFields: outputFields(C,F,P) vs. inputFields(C,F,P,D): dimensions C,F,P must match
189  const size_type f1[] = { 0, 1, 2 }, f2[] = { 0, 1, 2 };
190  for (size_type i=0; i<3; ++i) {
191  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
192  ">>> ERROR (ArrayTools::crossProductDataField): outputFields dimension does not match with inputFields");
193  }
194  } else {
195  // and rank-3 inputFields: outputFields(C,F,P) vs. inputFields(F,P,D): dimensions F,P must match
196  const size_type f1[] = { 1, 2 }, f2[] = { 0, 1 };
197  for (size_type i=0; i<2; ++i) {
198  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
199  ">>> ERROR (ArrayTools::crossProductDataField): outputFields dimension does not match with inputFields");
200  }
201  }
202  } else {
203  // In 3D:
204  if (inputFields.rank() == 4) {
205  // and rank-4 inputFields: outputFields(C,F,P,D) vs. inputFields(C,F,P,D): all dimensions C,F,P,D must match
206  for (size_type i=0; i<outputFields.rank(); ++i) {
207  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(i) != inputFields.extent(i), std::invalid_argument,
208  ">>> ERROR (ArrayTools::crossProductDataField): outputFields dimension does not match with inputFields");
209  }
210  } else {
211  // and rank-3 inputFields: outputFields(C,F,P,D) vs. inputFields(F,P,D): dimensions F,P,D must match
212  const size_type f1[] = { 1, 2, 3 }, f2[] = { 0, 1, 2 };
213  for (size_type i=0; i<3; ++i) {
214  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
215  ">>> ERROR (ArrayTools::crossProductDataField): outputFields dimension does not match with inputFields");
216  }
217  }
218  }
219  }
220 #endif
221  // body
222  ArrayTools<DeviceType>::Internal::crossProduct( outputFields,
223  inputData,
224  inputFields,
225  true );
226  }
227 
228 
229  template<typename DeviceType>
230  template<typename outputDataValueType, class ...outputDataProperties,
231  typename inputDataLeftValueType, class ...inputDataLeftProperties,
232  typename inputDataRightValueType, class ...inputDataRightProperties>
233  void
234  ArrayTools<DeviceType>::
235  crossProductDataData( Kokkos::DynRankView<outputDataValueType, outputDataProperties...> outputData,
236  const Kokkos::DynRankView<inputDataLeftValueType, inputDataLeftProperties...> inputDataLeft,
237  const Kokkos::DynRankView<inputDataRightValueType,inputDataRightProperties...> inputDataRight ) {
238 
239 #ifdef HAVE_INTREPID2_DEBUG
240  {
241  /*
242  * Check array rank and spatial dimension range (if applicable)
243  * (1) inputDataLeft(C,P,D);
244  * (2) inputDataRight(C,P,D) or (P,D);
245  * (3) outputData(C,P,D) in 3D, or (C,P) in 2D
246  */
247  // (1) inputDataLeft is (C, P, D) and 2 <= D <= 3 is required
248  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.rank() != 3, std::invalid_argument,
249  ">>> ERROR (ArrayTools::crossProductDataData): inputDataLeft must have rank 3");
250  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(2) < 2 || inputDataLeft.extent(2) > 3, std::invalid_argument,
251  ">>> ERROR (ArrayTools::crossProductDataData): inputDataLeft dimension(2) must be 2 or 3");
252 
253  // (2) inputDataRight is (C, P, D) or (P, D) and 2 <= (D=dimension(rank - 1)) <= 3 is required.
254  INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.rank() < 2 || inputDataRight.rank() > 3, std::invalid_argument,
255  ">>> ERROR (ArrayTools::crossProductDataData): inputDataRight must have rank 2 or 3" );
256  INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.extent(inputDataRight.rank()-1) < 2 ||
257  inputDataRight.extent(inputDataRight.rank()-1) > 3, std::invalid_argument,
258  ">>> ERROR (ArrayTools::crossProductDataData): inputDataRight dimension (rank-1) must have rank 2 or 3" );
259 
260  // (3) outputData is (C,P,D) in 3D and (C,P) in 2D => rank = inputDataLeft.extent(2)
261  INTREPID2_TEST_FOR_EXCEPTION( outputData.rank() != inputDataLeft.extent(2), std::invalid_argument,
262  ">>> ERROR (ArrayTools::crossProductDataData): outputData rank must match to inputDataLeft dimension(2)");
263  /*
264  * Dimension cross-checks:
265  * (1) inputDataLeft vs. inputDataRight
266  * (2) outputData vs. inputDataLeft
267  * (3) outputData vs. inputDataRight
268  *
269  * Cross-check (1):
270  */
271  if (inputDataRight.rank() == 3) {
272  // inputDataLeft(C,P,D) vs. inputDataRight(C,P,D): all dimensions C, P, D must match
273  for (size_type i=0; i<inputDataLeft.rank(); ++i) {
274  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(i) != inputDataRight.extent(i), std::invalid_argument,
275  ">>> ERROR (ArrayTools::crossProductDataData): inputDataLeft dimension does not match to inputDataRight");
276  }
277  }
278  // inputDataLeft(C, P,D) vs. inputDataRight(P,D): dimensions P, D must match
279  else {
280  const size_type f1[] = { 1, 2 }, f2[] = { 0, 1 };
281  for (size_type i=0; i<2; ++i) {
282  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
283  ">>> ERROR (ArrayTools::crossProductDataData): inputDataLeft dimension does not match to inputDataRight");
284  }
285  }
286  /*
287  * Cross-check (2):
288  */
289  if (inputDataLeft.extent(2) == 2) {
290  // in 2D: outputData(C,P) vs. inputDataLeft(C,P,D): dimensions C, P must match
291  const size_type f1[] = { 1, 2 }, f2[] = { 0, 1 };
292  for (size_type i=0; i<2; ++i) {
293  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != outputData.extent(f2[i]), std::invalid_argument,
294  ">>> ERROR (ArrayTools::crossProductDataData): inputDataLeft dimension does not match to outputData");
295  }
296  } else {
297  // in 3D: outputData(C,P,D) vs. inputDataLeft(C,P,D): all dimensions C, P, D must match
298  for (size_type i=0; i<inputDataLeft.rank(); ++i) {
299  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(i) != outputData.extent(i), std::invalid_argument,
300  ">>> ERROR (ArrayTools::crossProductDataData): inputDataLeft dimension does not match to outputData");
301  }
302  }
303  /*
304  * Cross-check (3):
305  */
306  if (inputDataLeft.extent(2) == 2) {
307  // In 2D:
308  if (inputDataRight.rank() == 3) {
309  // and rank-3 inputDataRight: outputData(C,P) vs. inputDataRight(C,P,D): dimensions C,P must match
310  const size_type f1[] = { 0, 1 }, f2[] = { 0, 1 };
311  for (size_type i=0; i<2; ++i) {
312  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
313  ">>> ERROR (ArrayTools::crossProductDataData): outputData dimension does not match to inputDataRight");
314  }
315  } else {
316  // and rank-2 inputDataRight: outputData(C,P) vs. inputDataRight(P,D): dimension P must match
317  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(1) != inputDataRight.extent(0), std::invalid_argument,
318  ">>> ERROR (ArrayTools::crossProductDataData): outputData dimension does not match to inputDataRight");
319  }
320  } else {
321  // In 3D:
322  if (inputDataRight.rank() == 3) {
323  // and rank-3 inputDataRight: outputData(C,P,D) vs. inputDataRight(C,P,D): all dimensions C,P,D must match
324  for (size_type i=0; i<outputData.rank(); ++i) {
325  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(i) != inputDataRight.extent(i), std::invalid_argument,
326  ">>> ERROR (ArrayTools::crossProductDataData): outputData dimension does not match to inputDataRight");
327  }
328  } else {
329  // and rank-2 inputDataRight: outputData(C,P,D) vs. inputDataRight(P,D): dimensions P, D must match
330  const size_type f1[] = { 1, 2 }, f2[] = { 0, 1 };
331  for (size_type i=0; i<2; ++i) {
332  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
333  ">>> ERROR (ArrayTools::crossProductDataData): outputData dimension does not match to inputDataRight");
334  }
335  }
336  }
337  }
338 #endif
339  // body
340  ArrayTools<DeviceType>::Internal::crossProduct( outputData,
341  inputDataLeft,
342  inputDataRight,
343  false );
344  }
345 
346 
347  namespace FunctorArrayTools {
351  template < typename OutputViewType, typename leftInputViewType, typename rightInputViewType >
352  struct F_outerProduct {
353  OutputViewType _output;
354  const leftInputViewType _leftInput;
355  const rightInputViewType _rightInput;
356  const bool _hasField;
357 
358  KOKKOS_INLINE_FUNCTION
359  F_outerProduct(OutputViewType output_,
360  leftInputViewType leftInput_,
361  rightInputViewType rightInput_,
362  const bool hasField_)
363  : _output(output_), _leftInput(leftInput_), _rightInput(rightInput_),
364  _hasField(hasField_) {}
365 
366  KOKKOS_INLINE_FUNCTION
367  void operator()(const size_type iter) const {
368  size_type cell, field(0), point;
369  size_type rightRank = _rightInput.rank();
370 
371  if (_hasField)
372  unrollIndex( cell, field, point,
373  _output.extent(0),
374  _output.extent(1),
375  _output.extent(2),
376  iter );
377  else
378  unrollIndex( cell, point,
379  _output.extent(0),
380  _output.extent(1),
381  iter );
382 
383  auto result = ( _hasField ? Kokkos::subview(_output, cell, field, point, Kokkos::ALL(), Kokkos::ALL()) :
384  Kokkos::subview(_output, cell, point, Kokkos::ALL(), Kokkos::ALL()));
385 
386  auto left = Kokkos::subview(_leftInput, cell, point, Kokkos::ALL());
387 
388  auto right = (rightRank == 2 + size_type(_hasField)) ?
389  ( _hasField ? Kokkos::subview(_rightInput, field, point, Kokkos::ALL()) :
390  Kokkos::subview(_rightInput, point, Kokkos::ALL())) :
391  ( _hasField ? Kokkos::subview(_rightInput, cell, field, point, Kokkos::ALL()) :
392  Kokkos::subview(_rightInput, cell, point, Kokkos::ALL()));
393 
394  const ordinal_type iend = result.extent(0);
395  const ordinal_type jend = result.extent(1);
396  for (ordinal_type i=0; i<iend; ++i)
397  for (ordinal_type j=0; j<jend; ++j)
398  result(i, j) = left(i)*right(j);
399  }
400  };
401  }
402 
403  template<typename DeviceType>
404  template<typename outputValueType, class ...outputProperties,
405  typename leftInputValueType, class ...leftInputProperties,
406  typename rightInputValueType, class ...rightInputProperties>
407  void
408  ArrayTools<DeviceType>::Internal::
409  outerProduct( Kokkos::DynRankView<outputValueType, outputProperties...> output,
410  const Kokkos::DynRankView<leftInputValueType, leftInputProperties...> leftInput,
411  const Kokkos::DynRankView<rightInputValueType,rightInputProperties...> rightInput,
412  const bool hasField ) {
413 
414  typedef Kokkos::DynRankView<outputValueType, outputProperties...> OutputViewType;
415  typedef const Kokkos::DynRankView<leftInputValueType, leftInputProperties...> leftInputViewType;
416  typedef const Kokkos::DynRankView<rightInputValueType,rightInputProperties...> rightInputViewType;
417  typedef FunctorArrayTools::F_outerProduct<OutputViewType, leftInputViewType, rightInputViewType> FunctorType;
418 
419  const size_type loopSize = ( hasField ? output.extent(0)*output.extent(1)*output.extent(2) :
420  output.extent(0)*output.extent(1) );
421  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
422  Kokkos::parallel_for( policy, FunctorType(output, leftInput, rightInput, hasField) );
423  }
424 
425 
426  template<typename DeviceType>
427  template<typename outputFieldValueType, class ...outputFieldProperties,
428  typename inputDataValueType, class ...inputDataProperties,
429  typename inputFieldValueType, class ...inputFieldProperties>
430  void
431  ArrayTools<DeviceType>::
432  outerProductDataField( Kokkos::DynRankView<outputFieldValueType,outputFieldProperties...> outputFields,
433  const Kokkos::DynRankView<inputDataValueType, inputDataProperties...> inputData,
434  const Kokkos::DynRankView<inputFieldValueType, inputFieldProperties...> inputFields ) {
435 
436 #ifdef HAVE_INTREPID2_DEBUG
437  {
438  /*
439  * Check array rank and spatial dimension range (if applicable)
440  * (1) inputData(C,P,D);
441  * (2) inputFields(C,F,P,D) or (F,P,D);
442  * (3) outputFields(C,F,P,D,D)
443  */
444  // (1) inputData is (C, P, D) and 2 <= D <= 3 is required
445  INTREPID2_TEST_FOR_EXCEPTION( inputData.rank() != 3, std::invalid_argument,
446  ">>> ERROR (ArrayTools::outerProductDataField): inputData must have rank 3");
447  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(2) < 2 || inputData.extent(2) > 3, std::invalid_argument,
448  ">>> ERROR (ArrayTools::outerProductDataField): inputData dimension(2) must be 2 or 3");
449 
450  // (2) inputFields is (C, F, P, D) or (F, P, D) and 2 <= (D=dimension(rank - 1)) <= 3 is required.
451  INTREPID2_TEST_FOR_EXCEPTION( inputFields.rank() < 3 || inputFields.rank() > 4, std::invalid_argument,
452  ">>> ERROR (ArrayTools::outerProductDataField): inputFields must have rank 3 or 4" );
453  INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(inputFields.rank()-1) < 2 ||
454  inputFields.extent(inputFields.rank()-1) < 3, std::invalid_argument,
455  ">>> ERROR (ArrayTools::outerProductDataField): inputFields dimension (rank-1) must have rank 2 or 3" );
456 
457  // (3) outputFields is (C,F,P,D,D)
458  INTREPID2_TEST_FOR_EXCEPTION( outputFields.rank() != 5, std::invalid_argument,
459  ">>> ERROR (ArrayTools::outerProductDataField): outputFields must have rank 5");
460  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(3) < 2 ||
461  outputFields.extent(3) > 3, std::invalid_argument,
462  ">>> ERROR (ArrayTools::outerProductDataField): outputFields dimension(3) must be 2 or 3");
463  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(4) < 2 ||
464  outputFields.extent(4) > 3, std::invalid_argument,
465  ">>> ERROR (ArrayTools::outerProductDataField): outputFields dimension(4) must be 2 or 3");
466 
467  /*
468  * Dimension cross-checks:
469  * (1) inputData vs. inputFields
470  * (2) outputFields vs. inputData
471  * (3) outputFields vs. inputFields
472  *
473  * Cross-check (2): outputFields(C,F,P,D,D) vs. inputData(C,P,D): dimensions C, P, D must match
474  */
475  {
476  const size_type f1[] = { 0, 2, 3, 4 }, f2[] = { 0, 1, 2, 2 };
477  for (size_type i=0; i<4; ++i) {
478  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputData.extent(f2[i]), std::invalid_argument,
479  ">>> ERROR (ArrayTools::outerProductDataField): outputFields dimension does not match with inputData");
480  }
481  }
482 
483  /*
484  * Cross-checks (1,3):
485  */
486  if (inputFields.rank() == 4) {
487  // Cross-check (1): inputData(C,P,D) vs. inputFields(C,F,P,D): dimensions C, P, D must match
488  {
489  const size_type f1[] = { 0, 1, 2 }, f2[] = { 0, 2, 3 };
490  for (size_type i=0; i<3; ++i) {
491  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
492  ">>> ERROR (ArrayTools::outerProductDataField): inputData dimension does not match with inputFields");
493  }
494  }
495 
496  // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(C,F,P,D): dimensions C, F, P, D must match
497  {
498  const size_type f1[] = { 0, 1, 2, 3, 4 }, f2[] = { 0, 1, 2, 3, 3 };
499  for (size_type i=0; i<5; ++i) {
500  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
501  ">>> ERROR (ArrayTools::outerProductDataField): outputFields dimension does not match with inputFields");
502  }
503  }
504  } else {
505  // Cross-check (1): inputData(C,P,D) vs. inputFields(F,P,D): dimensions P, D must match
506  {
507  const size_type f1[] = { 1, 2 }, f2[] = { 1, 2 };
508  for (size_type i=0; i<2; ++i) {
509  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
510  ">>> ERROR (ArrayTools::outerProductDataField): inputData dimension does not match with inputFields");
511  }
512  }
513 
514  // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(F,P,D): dimensions F, P, D must match
515  {
516  const size_type f1[] = { 1, 2, 3, 4 }, f2[] = { 0, 1, 2, 2 };
517  for (size_type i=0; i<4; ++i) {
518  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
519  ">>> ERROR (ArrayTools::outerProductDataField): outputFields dimension does not match with inputFields");
520  }
521  }
522  }
523  }
524 #endif
525  // body
526  ArrayTools<DeviceType>::Internal::outerProduct( outputFields,
527  inputData,
528  inputFields,
529  true );
530  }
531 
532 
533  template<typename DeviceType>
534  template<typename outputDataValueType, class ...outputDataProperties,
535  typename inputDataLeftValuetype, class ...inputDataLeftProperties,
536  typename inputDataRightValueType, class ...inputDataRightProperties>
537  void
538  ArrayTools<DeviceType>::
539  outerProductDataData( Kokkos::DynRankView<outputDataValueType, outputDataProperties...> outputData,
540  const Kokkos::DynRankView<inputDataLeftValuetype, inputDataLeftProperties...> inputDataLeft,
541  const Kokkos::DynRankView<inputDataRightValueType,inputDataRightProperties...> inputDataRight ) {
542 
543 #ifdef HAVE_INTREPID2_DEBUG
544  {
545  /*
546  * Check array rank and spatial dimension range (if applicable)
547  * (1) inputDataLeft(C,P,D);
548  * (2) inputDataRight(C,P,D) or (P,D);
549  * (3) outputData(C,P,D,D)
550  */
551  // (1) inputDataLeft is (C, P, D) and 2 <= D <= 3 is required
552  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.rank() != 3, std::invalid_argument,
553  ">>> ERROR (ArrayTools::outerProductDataField): inputDataLeft must have rank 3");
554  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(2) < 2 || inputDataLeft.extent(2) > 3, std::invalid_argument,
555  ">>> ERROR (ArrayTools::outerProductDataField): inputDataLeft dimension(2) must be 2 or 3");
556 
557  // (2) inputDataRight is (C, P, D) or (P, D) and 2 <= (D=dimension(rank - 1)) <= 3 is required.
558  INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.rank() < 2 || inputDataRight.rank() > 3, std::invalid_argument,
559  ">>> ERROR (ArrayTools::outerProductDataField): inputDataRight must have rank 2 or 3" );
560  INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.extent(inputDataRight.rank()-1) < 2 ||
561  inputDataRight.extent(inputDataRight.rank()-1) < 3, std::invalid_argument,
562  ">>> ERROR (ArrayTools::outerProductDataField): inputDataRight dimension (rank-1) must have rank 2 or 3" );
563 
564  // (3) outputData is (C,P,D,D)
565  INTREPID2_TEST_FOR_EXCEPTION( outputData.rank() != 4, std::invalid_argument,
566  ">>> ERROR (ArrayTools::outerProductDataField): outputData must have rank 5");
567  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(2) < 2 ||
568  outputData.extent(2) > 3, std::invalid_argument,
569  ">>> ERROR (ArrayTools::outerProductDataField): outputData dimension(3) must be 2 or 3");
570  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(3) < 2 ||
571  outputData.extent(3) > 3, std::invalid_argument,
572  ">>> ERROR (ArrayTools::outerProductDataField): outputData dimension(4) must be 2 or 3");
573 
574  /*
575  * Dimension cross-checks:
576  * (1) inputDataLeft vs. inputDataRight
577  * (2) outputData vs. inputDataLeft
578  * (3) outputData vs. inputDataRight
579  *
580  * Cross-check (2): outputData(C,P,D,D) vs. inputDataLeft(C,P,D): dimensions C, P, D must match
581  */
582  {
583  const size_type f1[] = { 0, 1, 2, 3 }, f2[] = { 0, 1, 2, 2 };
584  for (size_type i=0; i<4; ++i) {
585  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataLeft.extent(f2[i]), std::invalid_argument,
586  ">>> ERROR (ArrayTools::outerProductDataField): outputData dimension does not match with inputDataLeft");
587  }
588  }
589 
590  /*
591  * Cross-checks (1,3):
592  */
593  if (inputDataRight.rank() == 3) {
594  // Cross-check (1): inputDataLeft(C,P,D) vs. inputDataRight(C,P,D): all dimensions C, P, D must match
595  for (size_type i=0; i<inputDataLeft.rank(); ++i) {
596  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(i) != inputDataRight.extent(i), std::invalid_argument,
597  ">>> ERROR (ArrayTools::outerProductDataField): inputDataLeft dimension does not match with inputDataRight");
598  }
599 
600  // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(C,P,D): dimensions C, P, D must match
601  {
602  const size_type f1[] = { 0, 1, 2, 3 }, f2[] = { 0, 1, 2, 2 };
603  for (size_type i=0; i<4; ++i) {
604  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
605  ">>> ERROR (ArrayTools::outerProductDataField): outputData dimension does not match with inputDataRight");
606  }
607  }
608  } else {
609  // Cross-check (1): inputDataLeft(C,P,D) vs. inputDataRight(P,D): dimensions P, D must match
610  {
611  const size_type f1[] = { 1, 2 }, f2[] = { 0, 1 };
612  for (size_type i=0; i<2; ++i) {
613  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
614  ">>> ERROR (ArrayTools::outerProductDataField): inputDataLeft dimension does not match with inputDataRight");
615  }
616  }
617 
618  // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(P,D): dimensions P, D must match
619  {
620  const size_type f1[] = { 1, 2, 3 }, f2[] = { 0, 1, 1 };
621  for (size_type i=0; i<3; ++i) {
622  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
623  ">>> ERROR (ArrayTools::outerProductDataField): outputData dimension does not match with inputDataRight");
624  }
625  }
626  }
627  }
628 #endif
629  // body
630  ArrayTools<DeviceType>::Internal::outerProduct( outputData,
631  inputDataLeft,
632  inputDataRight,
633  false );
634  }
635 
636 
637  namespace FunctorArrayTools {
641  template < typename OutputViewType,
642  typename leftInputViewType,
643  typename rightInputViewType,
644  ordinal_type leftInputRank,
645  ordinal_type rightInputRank,
646  bool hasField,
647  bool isTranspose>
648  struct F_matvecProduct {
649  OutputViewType _output;
650  const leftInputViewType _leftInput;
651  const rightInputViewType _rightInput;
652 
653  const ordinal_type _iend;
654  const ordinal_type _jend;
655 
656  using value_type = typename OutputViewType::value_type;
657 
658  KOKKOS_INLINE_FUNCTION
659  F_matvecProduct(OutputViewType output_,
660  leftInputViewType leftInput_,
661  rightInputViewType rightInput_)
662  : _output(output_), _leftInput(leftInput_), _rightInput(rightInput_),
663  _iend(output_.extent_int(rank(output_)-1)), _jend(rightInput_.extent_int(rightInputRank-1))
664  {}
665 
666  // ****** hasField == true cases ******
667  KOKKOS_INLINE_FUNCTION
668  void operator()(const ordinal_type cl,
669  const ordinal_type bf,
670  const ordinal_type pt) const
671  {
672  apply_matvec_product(cl, bf, pt);
673  }
674 
675  template <ordinal_type l=leftInputRank, ordinal_type r=rightInputRank, bool hf=hasField>
676  KOKKOS_FORCEINLINE_FUNCTION
677  typename std::enable_if_t<l==4 && r==4 && hf, void>
678  apply_matvec_product(const ordinal_type &cl,
679  const ordinal_type &bf,
680  const ordinal_type &pt) const {
681  const auto lpt = (_leftInput.extent(1) == 1 ? size_type(0) : pt);
682  if (isTranspose) {
683  for (ordinal_type i=0;i<_iend;++i) {
684  value_type tmp(0);
685  for (ordinal_type j=0;j<_jend;++j)
686  tmp += _leftInput(cl,lpt, j,i)*_rightInput(cl,bf,pt, j);
687  _output(cl,bf,pt, i) = tmp;
688  }
689  } else {
690  for (ordinal_type i=0;i<_iend;++i) {
691  value_type tmp(0);
692  for (ordinal_type j=0;j<_jend;++j)
693  tmp += _leftInput(cl,lpt, i,j)*_rightInput(cl,bf,pt, j);
694  _output(cl,bf,pt, i) = tmp;
695  }
696  }
697  }
698 
699  template <ordinal_type l=leftInputRank, ordinal_type r=rightInputRank, bool hf=hasField>
700  KOKKOS_FORCEINLINE_FUNCTION
701  typename std::enable_if_t<l==4 && r==3 && hf, void>
702  apply_matvec_product(const ordinal_type &cl,
703  const ordinal_type &bf,
704  const ordinal_type &pt) const {
705  const auto lpt = (_leftInput.extent(1) == 1 ? size_type(0) : pt);
706  if (isTranspose) {
707  for (ordinal_type i=0;i<_iend;++i) {
708  value_type tmp(0);
709  for (ordinal_type j=0;j<_jend;++j)
710  tmp += _leftInput(cl,lpt, j,i)*_rightInput(bf,pt, j);
711  _output(cl,bf,pt, i) = tmp;
712  }
713  } else {
714  for (ordinal_type i=0;i<_iend;++i) {
715  value_type tmp(0);
716  for (ordinal_type j=0;j<_jend;++j)
717  tmp += _leftInput(cl,lpt, i,j)*_rightInput(bf,pt, j);
718  _output(cl,bf,pt, i) = tmp;
719  }
720  }
721  }
722 
723  template <ordinal_type l=leftInputRank, ordinal_type r=rightInputRank, bool hf=hasField>
724  KOKKOS_FORCEINLINE_FUNCTION
725  typename std::enable_if_t<l==3 && r==4 && hf, void>
726  apply_matvec_product(const ordinal_type &cl,
727  const ordinal_type &bf,
728  const ordinal_type &pt) const {
729  const auto lpt = (_leftInput.extent(1) == 1 ? size_type(0) : pt);
730  for (ordinal_type i=0;i<_iend;++i)
731  _output(cl,bf,pt, i) = _leftInput(cl,lpt, i)*_rightInput(cl,bf,pt, i);
732  }
733 
734  template <ordinal_type l=leftInputRank, ordinal_type r=rightInputRank, bool hf=hasField>
735  KOKKOS_FORCEINLINE_FUNCTION
736  typename std::enable_if_t<l==3 && r==3 && hf, void>
737  apply_matvec_product(const ordinal_type &cl,
738  const ordinal_type &bf,
739  const ordinal_type &pt) const {
740  const auto lpt = (_leftInput.extent(1) == 1 ? size_type(0) : pt);
741  for (ordinal_type i=0;i<_iend;++i)
742  _output(cl,bf,pt, i) = _leftInput(cl,lpt, i)*_rightInput(bf,pt, i);
743  }
744 
745  template <ordinal_type l=leftInputRank, ordinal_type r=rightInputRank, bool hf=hasField>
746  KOKKOS_FORCEINLINE_FUNCTION
747  typename std::enable_if_t<l==2 && r==4 && hf, void>
748  apply_matvec_product(const ordinal_type &cl,
749  const ordinal_type &bf,
750  const ordinal_type &pt) const {
751  const auto lpt = (_leftInput.extent(1) == 1 ? size_type(0) : pt);
752  const value_type & val = _leftInput(cl,lpt);
753  for (ordinal_type i=0;i<_iend;++i) {
754  _output(cl,bf,pt, i) = val*_rightInput(cl,bf,pt, i);
755  }
756  }
757 
758  template <ordinal_type l=leftInputRank, ordinal_type r=rightInputRank, bool hf=hasField>
759  KOKKOS_FORCEINLINE_FUNCTION
760  typename std::enable_if_t<l==2 && r==3 && hf, void>
761  apply_matvec_product(const ordinal_type &cl,
762  const ordinal_type &bf,
763  const ordinal_type &pt) const {
764  const auto lpt = (_leftInput.extent(1) == 1 ? size_type(0) : pt);
765  const value_type & val = _leftInput(cl,lpt);
766  for (ordinal_type i=0;i<_iend;++i) {
767  _output(cl,bf,pt, i) = val*_rightInput(bf,pt, i);
768  }
769  }
770 
771  // ****** hasField == false cases ******
772  KOKKOS_INLINE_FUNCTION
773  void operator()(const ordinal_type cl,
774  const ordinal_type pt) const
775  {
776  apply_matvec_product(cl, pt);
777  }
778 
779  template <ordinal_type l=leftInputRank, ordinal_type r=rightInputRank, bool hf=hasField>
780  KOKKOS_FORCEINLINE_FUNCTION
781  typename std::enable_if_t<l==4 && r==3 && !hf, void>
782  apply_matvec_product(const ordinal_type &cl,
783  const ordinal_type &pt) const {
784  const auto lpt = (_leftInput.extent(1) == 1 ? size_type(0) : pt);
785  if (isTranspose) {
786  for (ordinal_type i=0;i<_iend;++i) {
787  value_type tmp(0);
788  for (ordinal_type j=0;j<_jend;++j)
789  tmp += _leftInput(cl,lpt, j,i)*_rightInput(cl,pt, j);
790  _output(cl,pt, i) = tmp;
791  }
792  } else {
793  for (ordinal_type i=0;i<_iend;++i) {
794  value_type tmp(0);
795  for (ordinal_type j=0;j<_jend;++j)
796  tmp += _leftInput(cl,lpt, i,j)*_rightInput(cl,pt, j);
797  _output(cl,pt, i) = tmp;
798  }
799  }
800  }
801 
802  template <ordinal_type l=leftInputRank, ordinal_type r=rightInputRank, bool hf=hasField>
803  KOKKOS_FORCEINLINE_FUNCTION
804  typename std::enable_if_t<l==4 && r==2 && !hf, void>
805  apply_matvec_product(const ordinal_type &cl,
806  const ordinal_type &pt) const {
807  const auto lpt = (_leftInput.extent(1) == 1 ? size_type(0) : pt);
808  if (isTranspose) {
809  for (ordinal_type i=0;i<_iend;++i) {
810  value_type tmp(0);
811  for (ordinal_type j=0;j<_jend;++j)
812  tmp += _leftInput(cl,lpt, j,i)*_rightInput(pt, j);
813  _output(cl,pt, i) = tmp;
814  }
815  } else {
816  for (ordinal_type i=0;i<_iend;++i) {
817  value_type tmp(0);
818  for (ordinal_type j=0;j<_jend;++j)
819  tmp += _leftInput(cl,lpt, i,j)*_rightInput(pt, j);
820  _output(cl,pt, i) = tmp;
821  }
822  }
823  }
824 
825  template <ordinal_type l=leftInputRank, ordinal_type r=rightInputRank, bool hf=hasField>
826  KOKKOS_FORCEINLINE_FUNCTION
827  typename std::enable_if_t<l==3 && r==3 && !hf, void>
828  apply_matvec_product(const ordinal_type &cl,
829  const ordinal_type &pt) const {
830  const auto lpt = (_leftInput.extent(1) == 1 ? size_type(0) : pt);
831  for (ordinal_type i=0;i<_iend;++i)
832  _output(cl,pt, i) = _leftInput(cl,lpt, i)*_rightInput(cl,pt, i);
833  }
834 
835  template <ordinal_type l=leftInputRank, ordinal_type r=rightInputRank, bool hf=hasField>
836  KOKKOS_FORCEINLINE_FUNCTION
837  typename std::enable_if_t<l==3 && r==2 && !hf, void>
838  apply_matvec_product(const ordinal_type &cl,
839  const ordinal_type &pt) const {
840  const auto lpt = (_leftInput.extent(1) == 1 ? size_type(0) : pt);
841  for (ordinal_type i=0;i<_iend;++i)
842  _output(cl,pt, i) = _leftInput(cl,lpt, i)*_rightInput(pt, i);
843  }
844 
845  template <ordinal_type l=leftInputRank, ordinal_type r=rightInputRank, bool hf=hasField>
846  KOKKOS_FORCEINLINE_FUNCTION
847  typename std::enable_if_t<l==2 && r==3 && !hf, void>
848  apply_matvec_product(const ordinal_type &cl,
849  const ordinal_type &pt) const {
850  const auto lpt = (_leftInput.extent(1) == 1 ? size_type(0) : pt);
851  const value_type & val = _leftInput(cl,lpt);
852  for (ordinal_type i=0;i<_iend;++i) {
853  _output(cl,pt, i) = val*_rightInput(cl,pt, i);
854  }
855  }
856 
857  template <ordinal_type l=leftInputRank, ordinal_type r=rightInputRank, bool hf=hasField>
858  KOKKOS_FORCEINLINE_FUNCTION
859  typename std::enable_if_t<l==2 && r==2 && !hf, void>
860  apply_matvec_product(const ordinal_type &cl,
861  const ordinal_type &pt) const {
862  const auto lpt = (_leftInput.extent(1) == 1 ? size_type(0) : pt);
863  const value_type & val = _leftInput(cl,lpt);
864  for (ordinal_type i=0;i<_iend;++i) {
865  _output(cl,pt, i) = val*_rightInput(pt, i);
866  }
867  }
868  };
869  } //namespace
870 
871  template<typename DeviceType>
872  template<typename outputValueType, class ...outputProperties,
873  typename leftInputValueType, class ...leftInputProperties,
874  typename rightInputValueType, class ...rightInputProperties>
875  void
876  ArrayTools<DeviceType>::Internal::
877  matvecProduct( Kokkos::DynRankView<outputValueType, outputProperties...> output,
878  const Kokkos::DynRankView<leftInputValueType, leftInputProperties...> leftInput,
879  const Kokkos::DynRankView<rightInputValueType,rightInputProperties...> rightInput,
880  const bool hasField,
881  const bool isTranspose ) {
882 
883  using Output = Kokkos::DynRankView<outputValueType, outputProperties...>;
884  using Left = const Kokkos::DynRankView<leftInputValueType, leftInputProperties...>;
885  using Right = const Kokkos::DynRankView<rightInputValueType,rightInputProperties...>;
886 
887  // FTNMAB: FunctorType with left rank N, right rank M, hasField = (A==T), isTranspose = (B==T)
888 
889  // hasField == true
890  using FT44TT = FunctorArrayTools::F_matvecProduct<Output,Left,Right,4,4, true, true>;
891  using FT43TT = FunctorArrayTools::F_matvecProduct<Output,Left,Right,4,3, true, true>;
892  using FT44TF = FunctorArrayTools::F_matvecProduct<Output,Left,Right,4,4, true, false>;
893  using FT43TF = FunctorArrayTools::F_matvecProduct<Output,Left,Right,4,3, true, false>;
894 
895  // for left rank 3, 2, isTranspose does not matter
896  using FT34TF = FunctorArrayTools::F_matvecProduct<Output,Left,Right,3,4, true, false>;
897  using FT33TF = FunctorArrayTools::F_matvecProduct<Output,Left,Right,3,3, true, false>;
898  using FT24TF = FunctorArrayTools::F_matvecProduct<Output,Left,Right,2,4, true, false>;
899  using FT23TF = FunctorArrayTools::F_matvecProduct<Output,Left,Right,2,3, true, false>;
900 
901  // hasField == false
902  using FT43FT = FunctorArrayTools::F_matvecProduct<Output,Left,Right,4,3, false, true>;
903  using FT42FT = FunctorArrayTools::F_matvecProduct<Output,Left,Right,4,2, false, true>;
904  using FT43FF = FunctorArrayTools::F_matvecProduct<Output,Left,Right,4,3, false, false>;
905  using FT42FF = FunctorArrayTools::F_matvecProduct<Output,Left,Right,4,2, false, false>;
906 
907  // for left rank 3, 2, isTranspose does not matter
908  using FT33FF = FunctorArrayTools::F_matvecProduct<Output,Left,Right,3,3, false, false>;
909  using FT32FF = FunctorArrayTools::F_matvecProduct<Output,Left,Right,3,2, false, false>;
910  using FT23FF = FunctorArrayTools::F_matvecProduct<Output,Left,Right,2,3, false, false>;
911  using FT22FF = FunctorArrayTools::F_matvecProduct<Output,Left,Right,2,2, false, false>;
912 
913  const ordinal_type l = leftInput.rank();
914  const ordinal_type r = rightInput.rank();
915 
916  using range_policy2 = Kokkos::MDRangePolicy< ExecSpaceType, Kokkos::Rank<2>, Kokkos::IndexType<ordinal_type> >;
917  range_policy2 policy2( { 0, 0 }, { output.extent(0), output.extent(1) } );
918 
919  using range_policy3 = Kokkos::MDRangePolicy< ExecSpaceType, Kokkos::Rank<3>, Kokkos::IndexType<ordinal_type> >;
920  range_policy3 policy3( { 0, 0, 0 }, { output.extent(0), output.extent(1), output.extent(2) } );
921 
922  // just to make the below a little easier to read, we pack l and r together:
923  const ordinal_type lr = l * 10 + r;
924  const auto &ov = output;
925  const auto &lv = leftInput;
926  const auto &rv = rightInput;
927 
928  if (hasField) // this means we want policy3
929  {
930  if (l == 4)
931  {
932  // isTranspose matters
933  if (isTranspose)
934  {
935  switch (r)
936  {
937  case 4: { FT44TT f(ov,lv,rv); Kokkos::parallel_for( policy3, f ); } break;
938  default: { FT43TT f(ov,lv,rv); Kokkos::parallel_for( policy3, f ); }
939  }
940  }
941  else
942  {
943  switch (r)
944  {
945  case 4: { FT44TF f(ov,lv,rv); Kokkos::parallel_for( policy3, f ); } break;
946  default: { FT43TF f(ov,lv,rv); Kokkos::parallel_for( policy3, f ); }
947  }
948  }
949  }
950  else // l == 3 or 2; isTranspose does not matter
951  {
952  switch (lr)
953  {
954  case 34: { FT34TF f(ov,lv,rv); Kokkos::parallel_for( policy3, f ); } break;
955  case 33: { FT33TF f(ov,lv,rv); Kokkos::parallel_for( policy3, f ); } break;
956  case 24: { FT24TF f(ov,lv,rv); Kokkos::parallel_for( policy3, f ); } break;
957  default: { FT23TF f(ov,lv,rv); Kokkos::parallel_for( policy3, f ); } // 23
958  }
959  }
960  }
961  else // hasField is false; use policy2
962  {
963  if (l == 4)
964  {
965  // isTranspose matters
966  if (isTranspose)
967  {
968  switch (r)
969  {
970  case 3: { FT43FT f(ov,lv,rv); Kokkos::parallel_for( policy2, f ); } break;
971  default: { FT42FT f(ov,lv,rv); Kokkos::parallel_for( policy2, f ); } // 2
972  }
973  }
974  else
975  {
976  switch (r)
977  {
978  case 3: { FT43FF f(ov,lv,rv); Kokkos::parallel_for( policy2, f ); } break;
979  default: { FT42FF f(ov,lv,rv); Kokkos::parallel_for( policy2, f ); } // 2
980  }
981  }
982  }
983  else // l == 3 or 2; isTranspose does not matter
984  {
985  switch (lr)
986  {
987  case 33: { FT33FF f(ov,lv,rv); Kokkos::parallel_for( policy2, f ); } break;
988  case 32: { FT32FF f(ov,lv,rv); Kokkos::parallel_for( policy2, f ); } break;
989  case 23: { FT23FF f(ov,lv,rv); Kokkos::parallel_for( policy2, f ); } break;
990  default: { FT22FF f(ov,lv,rv); Kokkos::parallel_for( policy2, f ); } // 22
991  }
992  }
993  }
994  }
995 
996  template<typename DeviceType>
997  template<typename outputFieldValueType, class ...outputFieldProperties,
998  typename inputDataValueType, class ...inputDataProperties,
999  typename inputFieldValueType, class ...inputFieldProperties>
1000  void
1001  ArrayTools<DeviceType>::matvecProductDataField( Kokkos::DynRankView<outputFieldValueType,outputFieldProperties...> outputFields,
1002  const Kokkos::DynRankView<inputDataValueType, inputDataProperties...> inputData,
1003  const Kokkos::DynRankView<inputFieldValueType, inputFieldProperties...> inputFields,
1004  const char transpose ) {
1005 
1006 #ifdef HAVE_INTREPID2_DEBUG
1007  {
1008  /*
1009  * Check array rank and spatial dimension range (if applicable)
1010  * (1) inputData(C,P), (C,P,D) or (C,P,D,D); P=1 is admissible to allow multiply by const. data
1011  * (2) inputFields(C,F,P,D) or (F,P,D);
1012  * (3) outputFields(C,F,P,D)
1013  */
1014  // (1) inputData is (C,P), (C, P, D) or (C, P, D, D) and 1 <= D <= 3 is required
1015  INTREPID2_TEST_FOR_EXCEPTION( inputData.rank() < 2 || inputData.rank() > 4, std::invalid_argument,
1016  ">>> ERROR (ArrayTools::matvecProductDataField): inputData must have rank 2,3 or 4" );
1017  if (inputData.rank() > 2) {
1018  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(2) < 1 ||
1019  inputData.extent(2) > 3, std::invalid_argument,
1020  ">>> ERROR (ArrayTools::matvecProductDataField): inputData dimension(2) must be 1,2 or 3");
1021  }
1022  if (inputData.rank() == 4) {
1023  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(3) < 1 ||
1024  inputData.extent(3) > 3, std::invalid_argument,
1025  ">>> ERROR (ArrayTools::matvecProductDataField): inputData dimension(3) must be 1,2 or 3");
1026  }
1027 
1028  // (2) inputFields is (C, F, P, D) or (F, P, D) and 1 <= (D=dimension(rank - 1)) <= 3 is required.
1029  INTREPID2_TEST_FOR_EXCEPTION( inputFields.rank() < 3 ||
1030  inputFields.rank() > 4, std::invalid_argument,
1031  ">>> ERROR (ArrayTools::matvecProductDataField): inputFields must have rank 3 or 4" );
1032  INTREPID2_TEST_FOR_EXCEPTION( (inputFields.rank()-1) < 1 ||
1033  (inputFields.rank()-1) > 3, std::invalid_argument,
1034  ">>> ERROR (ArrayTools::matvecProductDataField): inputFields dimension(rank-1) must be 1,2, or 3" );
1035 
1036  // (3) outputFields is (C,F,P,D)
1037  INTREPID2_TEST_FOR_EXCEPTION( outputFields.rank() != 4, std::invalid_argument,
1038  ">>> ERROR (ArrayTools::matvecProductDataField): outputFields must have rank 4" );
1039  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(3) < 1 ||
1040  outputFields.extent(3) > 3, std::invalid_argument,
1041  ">>> ERROR (ArrayTools::matvecProductDataField): outputFields dimension(3) must be 1,2 or 3" );
1042 
1043  /*
1044  * Dimension cross-checks:
1045  * (1) inputData vs. inputFields
1046  * (2) outputFields vs. inputData
1047  * (3) outputFields vs. inputFields
1048  *
1049  * Cross-check (2): outputFields(C,F,P,D) vs. inputData(C,P), (C,P,D) or (C,P,D,D):
1050  * dimensions C, and D must match in all cases, dimension P must match only when non-constant
1051  * data is specified (P>1). Do not check P dimensions with constant data, i.e., when P=1 in
1052  * inputData(C,1,...)
1053  */
1054  if (inputData.extent(1) > 1) { // check P dimension if P>1 in inputData
1055  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(2) != inputData.extent(1), std::invalid_argument,
1056  ">>> ERROR (ArrayTools::matvecProductDataField): outputFields dimension(2) must match to inputData dimension(1)" );
1057  }
1058  if (inputData.rank() == 2) { // inputData(C,P) -> C match
1059  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(0) != inputData.extent(0), std::invalid_argument,
1060  ">>> ERROR (ArrayTools::matvecProductDataField): outputFields dimension(0) must match to inputData dimension(0)" );
1061  }
1062  if (inputData.rank() == 3) { // inputData(C,P,D) -> C, D match
1063  const size_type f1[] = { 0, 3 }, f2[] = { 0, 2 };
1064  for (size_type i=0; i<2; ++i) {
1065  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputData.extent(f2[i]), std::invalid_argument,
1066  ">>> ERROR (ArrayTools::matvecProductDataField): outputFields dimension does not match to inputData dimension" );
1067  }
1068  }
1069  if (inputData.rank() == 4) { // inputData(C,P,D,D) -> C, D, D match
1070  const size_type f1[] = { 0, 3, 3 }, f2[] = { 0, 2, 3 };
1071  for (size_type i=0; i<3; ++i) {
1072  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputData.extent(f2[i]), std::invalid_argument,
1073  ">>> ERROR (ArrayTools::matvecProductDataField): outputFields dimension does not match to inputData dimension" );
1074  }
1075  }
1076 
1077  /*
1078  * Cross-checks (1,3):
1079  */
1080  if (inputFields.rank() == 4) {
1081  // Cross-check (1): inputData(C,P), (C,P,D) or (C,P,D,D) vs. inputFields(C,F,P,D): dimensions C, P, D must match
1082  if (inputData.extent(1) > 1) { // check P dimension if P>1 in inputData
1083  INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(2) != inputData.extent(1), std::invalid_argument,
1084  ">>> ERROR (ArrayTools::matvecProductDataField): inputFields dimension (2) does not match to inputData dimension (1)" );
1085  }
1086  if (inputData.rank() == 2) { // inputData(C,P) -> C match
1087  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(0) != inputFields.extent(0), std::invalid_argument,
1088  ">>> ERROR (ArrayTools::matvecProductDataField): inputData dimension (0) does not match to inputFields dimension (0)" );
1089  }
1090  if (inputData.rank() == 3) { // inputData(C,P,D) -> C, D match
1091  const size_type f1[] = { 0, 2 }, f2[] = { 0, 3 };
1092  for (size_type i=0; i<2; ++i) {
1093  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
1094  ">>> ERROR (ArrayTools::matvecProductDataField): inputData dimension does not match to inputFields dimension" );
1095  }
1096  }
1097  if (inputData.rank() == 4) { // inputData(C,P,D,D) -> C, D, D match
1098  const size_type f1[] = { 0, 2, 3 }, f2[] = { 0, 3, 3 };
1099  for (size_type i=0; i<3; ++i) {
1100  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
1101  ">>> ERROR (ArrayTools::matvecProductDataField): inputData dimension does not match to inputFields dimension" );
1102  }
1103  }
1104 
1105  // Cross-check (3): outputFields(C,F,P,D) vs. inputFields(C,F,P,D): all dimensions C, F, P, D must match
1106  for (size_type i=0; i<outputFields.rank(); ++i) {
1107  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(i) != inputFields.extent(i), std::invalid_argument,
1108  ">>> ERROR (ArrayTools::matvecProductDataField): outputFields dimension does not match to inputFields dimension" );
1109  }
1110  } else {
1111  // Cross-check (1): inputData(C,P), (C,P,D) or (C,P,D,D) vs. inputFields(F,P,D): dimensions P, D must match
1112  if (inputData.extent(1) > 1) { // check P if P>1 in inputData
1113  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(1) != inputFields.extent(1), std::invalid_argument,
1114  ">>> ERROR (ArrayTools::matvecProductDataField): inputData dimension(1) does not match to inputFields dimension (1)" );
1115  }
1116  if (inputData.rank() == 3) {
1117  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(2) != inputFields.extent(2), std::invalid_argument,
1118  ">>> ERROR (ArrayTools::matvecProductDataField): inputData dimension(2) does not match to inputFields dimension (2)" );
1119  }
1120  if (inputData.rank() == 4) {
1121  const size_type f1[] = { 2, 3 }, f2[] = { 2, 2 };
1122  for (size_type i=0; i<2; ++i) {
1123  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
1124  ">>> ERROR (ArrayTools::matvecProductDataField): inputData dimension does not match to inputFields dimension" );
1125  }
1126  }
1127 
1128  // Cross-check (3): outputFields(C,F,P,D) vs. inputFields(F,P,D): dimensions F, P, D must match
1129  {
1130  const size_type f1[] = { 1, 2, 3 }, f2[] = { 0, 1, 2 };
1131  for (size_type i=0; i<3; ++i) {
1132  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
1133  ">>> ERROR (ArrayTools::matvecProductDataField): outputFields dimension does not match to inputFields dimension" );
1134  }
1135  }
1136  }
1137  }
1138 #endif
1139  // body
1140  ArrayTools<DeviceType>::Internal::matvecProduct( outputFields,
1141  inputData,
1142  inputFields,
1143  true,
1144  transpose == 't' || transpose == 'T' );
1145  }
1146 
1147 
1148 
1149  template<typename DeviceType>
1150  template<typename outputDataValueType, class ...outputDataProperties,
1151  typename inputDataLeftValueType, class ...inputDataLeftProperties,
1152  typename inputDataRightValueType, class ...inputDataRightProperties>
1153  void
1154  ArrayTools<DeviceType>::
1155  matvecProductDataData( Kokkos::DynRankView<outputDataValueType, outputDataProperties...> outputData,
1156  const Kokkos::DynRankView<inputDataLeftValueType, inputDataLeftProperties...> inputDataLeft,
1157  const Kokkos::DynRankView<inputDataRightValueType,inputDataRightProperties...> inputDataRight,
1158  const char transpose ) {
1159 
1160 #ifdef HAVE_INTREPID2_DEBUG
1161  {
1162  /*
1163  * Check array rank and spatial dimension range (if applicable)
1164  * (1) inputDataLeft(C,P), (C,P,D) or (C,P,D1,D2); P=1 is admissible to allow multiply by const. left data
1165  * (2) inputDataRight(C,P,D) or (P,D);
1166  * (3) outputData(C,P,D)
1167  */
1168  // (1) inputDataLeft is (C,P), (C,P,D) or (C,P,D1,D2) and 1 <= D,D1,D2 <= 3 is required
1169  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.rank() < 2 ||
1170  inputDataLeft.rank() > 4, std::invalid_argument,
1171  ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft must have rank 2,3 or 4" );
1172 
1173  if (inputDataLeft.rank() > 2) {
1174  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(2) < 1 ||
1175  inputDataLeft.extent(2) > 3, std::invalid_argument,
1176  ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension(2) must be 1, 2 or 3");
1177  }
1178  if (inputDataLeft.rank() == 4) {
1179  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(3) < 1 ||
1180  inputDataLeft.extent(3) > 3, std::invalid_argument,
1181  ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension(3) must be 1, 2 or 3");
1182  }
1183 
1184  // (2) inputDataRight is (C, P, D) or (P, D) and 1 <= (D=dimension(rank - 1)) <= 3 is required.
1185  INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.rank() < 2 ||
1186  inputDataRight.rank() > 3, std::invalid_argument,
1187  ">>> ERROR (ArrayTools::matvecProductDataData): inputDataRight must have rank 2 or 3" );
1188  INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.extent(inputDataRight.rank()-1) < 1 ||
1189  inputDataRight.extent(inputDataRight.rank()-1) > 3, std::invalid_argument,
1190  ">>> ERROR (ArrayTools::matvecProductDataData): inputDataRight dimension (rank-1) must be 1,2 or 3" );
1191 
1192  // (3) outputData is (C,P,D)
1193  INTREPID2_TEST_FOR_EXCEPTION( outputData.rank() != 3, std::invalid_argument,
1194  ">>> ERROR (ArrayTools::matvecProductDataData): outputData must have rank 3" );
1195  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(2) < 1 ||
1196  outputData.extent(2) > 3, std::invalid_argument,
1197  ">>> ERROR (ArrayTools::matvecProductDataData): outputData dimension(2) must be 1, 2 or 3");
1198 
1199  /*
1200  * Dimension cross-checks:
1201  * (1) inputDataLeft vs. inputDataRight
1202  * (2) outputData vs. inputDataLeft
1203  * (3) outputData vs. inputDataRight
1204  *
1205  * Cross-check (2): outputData(C,P,D) vs. inputDataLeft(C,P), (C,P,D) or (C,P,D,D):
1206  * dimensions C, and D must match in all cases, dimension P must match only when non-constant
1207  * data is specified (P>1). Do not check P dimensions with constant left data, i.e., when P=1 in
1208  * inputDataLeft(C,1,...)
1209  */
1210  if (inputDataLeft.extent(1) > 1) { // check P dimension if P>1 in inputDataLeft
1211  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(1) != inputDataLeft.extent(1), std::invalid_argument,
1212  ">>> ERROR (ArrayTools::matvecProductDataData): outputData dimension(1) muat match inputDataLeft dimension(1)" );
1213  }
1214  if (inputDataLeft.rank() == 2) { // inputDataLeft(C,P): check C
1215  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(0) != inputDataLeft.extent(0), std::invalid_argument,
1216  ">>> ERROR (ArrayTools::matvecProductDataData): outputData dimension(0) muat match inputDataLeft dimension(0)" );
1217  }
1218  if (inputDataLeft.rank() == 3) { // inputDataLeft(C,P,D): check C and D
1219  const size_type f1[] = { 0, 2 }, f2[] = { 0, 2 };
1220  for (size_type i=0; i<2; ++i) {
1221  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataLeft.extent(f2[i]), std::invalid_argument,
1222  ">>> ERROR (ArrayTools::matvecProductDataData): outputData dimension muat match inputDataLeft dimension" );
1223  }
1224  }
1225  if (inputDataLeft.rank() == 4) { // inputDataLeft(C,P,D1,D2): check C and D
1226  size_type f1[] = { 0, 2}, f2[] = { 0, 2};
1227  if (transpose) f2[1] = 3;
1228  for (size_type i=0; i<2; ++i) {
1229  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataLeft.extent(f2[i]), std::invalid_argument,
1230  ">>> ERROR (ArrayTools::matvecProductDataData): outputData dimension muat match inputDataLeft dimension" );
1231  }
1232  }
1233 
1234  /*
1235  * Cross-checks (1,3):
1236  */
1237  if (inputDataRight.rank() == 3) {
1238  // Cross-check (1): inputDataLeft(C,P), (C,P,D), or (C,P,D,D) vs. inputDataRight(C,P,D): dimensions C, P, D must match
1239  if (inputDataLeft.extent(1) > 1) { // check P dimension if P>1 in inputDataLeft
1240  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(1) != inputDataRight.extent(1), std::invalid_argument,
1241  ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension (1) muat match to inputDataRight dimension (1)" );
1242  }
1243  if (inputDataLeft.rank() == 2) { // inputDataLeft(C,P): check C
1244  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(0) != inputDataRight.extent(0), std::invalid_argument,
1245  ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension (0) muat match to inputDataRight dimension (0)" );
1246  }
1247  if (inputDataLeft.rank() == 3) { // inputDataLeft(C,P,D): check C and D
1248  const size_type f1[] = { 0, 2 }, f2[] = { 0, 2 };
1249  for (size_type i=0; i<2; ++i) {
1250  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1251  ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension muat match to inputDataRight dimension" );
1252  }
1253  }
1254  if (inputDataLeft.rank() == 4) { // inputDataLeft(C,P,D,D): check C and D
1255  size_type f1[] = { 0, 3}, f2[] = { 0, 2};
1256  if (transpose) f1[1] = 2;
1257  for (size_type i=0; i<2; ++i) {
1258  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1259  ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension muat match to inputDataRight dimension" );
1260  }
1261  }
1262 
1263  // Cross-check (3): outputData(C,P) vs. inputDataRight(C,P): all dimensions C, P must match
1264  for (size_type i=0; i<outputData.rank()-1; ++i) {
1265  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(i) != inputDataRight.extent(i), std::invalid_argument,
1266  ">>> ERROR (ArrayTools::matvecProductDataData): outputData dimension muat match to inputDataRight dimension" );
1267  }
1268  } else {
1269  // Cross-check (1): inputDataLeft(C,P), (C,P,D), or (C,P,D,D) vs. inputDataRight(P,D): dimensions P, D must match
1270  if (inputDataLeft.extent(1) > 1) { // check P if P>1 in inputData
1271  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(1) != inputDataRight.extent(0), std::invalid_argument,
1272  ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension (1) does mot match to inputDataright dimension (0)" );
1273  }
1274  if (inputDataLeft.rank() == 3) { // inputDataLeft(C,P,D): check D
1275  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(2) != inputDataRight.extent(1), std::invalid_argument,
1276  ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension (2) does mot match to inputDataright dimension (1)" );
1277  }
1278  if (inputDataLeft.rank() == 4) { // inputDataLeft(C,P,D,D): check D
1279  const size_type f1[] = { 2, 3 }, f2[] = { 1, 1 };
1280  for (size_type i=0; i<2; ++i) {
1281  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1282  ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension muat match to inputDataRight dimension" );
1283  }
1284  }
1285 
1286  // Cross-check (3): outputData(C,P,D) vs. inputDataRight(P,D): dimensions P, D must match
1287  {
1288  const size_type f1[] = { 1, 2 }, f2[] = { 0, 1 };
1289  for (size_type i=0; i<2; ++i) {
1290  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1291  ">>> ERROR (ArrayTools::matvecProductDataData): outputData dimension muat match to inputDataRight dimension" );
1292  }
1293  }
1294  }
1295  }
1296 #endif
1297  // body
1298  ArrayTools<DeviceType>::Internal::matvecProduct( outputData,
1299  inputDataLeft,
1300  inputDataRight,
1301  false,
1302  transpose == 't' || transpose == 'T' );
1303  }
1304 
1305 
1306  namespace FunctorArrayTools {
1310  template < typename OutputViewType, typename leftInputViewType, typename rightInputViewType >
1311  struct F_matmatProduct{
1312  OutputViewType _output;
1313  leftInputViewType _leftInput;
1314  rightInputViewType _rightInput;
1315  const bool _hasField, _isTranspose;
1316  typedef typename leftInputViewType::value_type value_type;
1317 
1318  KOKKOS_INLINE_FUNCTION
1319  F_matmatProduct(OutputViewType output_,
1320  leftInputViewType leftInput_,
1321  rightInputViewType rightInput_,
1322  const bool hasField_,
1323  const bool isTranspose_)
1324  : _output(output_), _leftInput(leftInput_), _rightInput(rightInput_),
1325  _hasField(hasField_), _isTranspose(isTranspose_) {}
1326 
1327  KOKKOS_INLINE_FUNCTION
1328  void operator()(const size_type iter) const {
1329  size_type cell(0), field(0), point(0);
1330  size_type leftRank = _leftInput.rank();
1331  size_type rightRank = _rightInput.rank();
1332 
1333  if (_hasField)
1334  unrollIndex( cell, field, point,
1335  _output.extent(0),
1336  _output.extent(1),
1337  _output.extent(2),
1338  iter );
1339  else
1340  unrollIndex( cell, point,
1341  _output.extent(0),
1342  _output.extent(1),
1343  iter );
1344 
1345  auto result = ( _hasField ? Kokkos::subview(_output, cell, field, point, Kokkos::ALL(), Kokkos::ALL()) :
1346  Kokkos::subview(_output, cell, point, Kokkos::ALL(), Kokkos::ALL()));
1347 
1348  const auto lpoint = (_leftInput.extent(1) == 1 ? 0 : point);
1349  auto left = ( leftRank == 4 ? Kokkos::subview(_leftInput, cell, lpoint, Kokkos::ALL(), Kokkos::ALL()) :
1350  leftRank == 3 ? Kokkos::subview(_leftInput, cell, lpoint, Kokkos::ALL()) :
1351  Kokkos::subview(_leftInput, cell, lpoint) );
1352 
1353  //temporary
1354  const bool hasCell = (_hasField ? rightRank == 5 : rightRank == 4);
1355  auto right = ( _hasField ? ( hasCell ? Kokkos::subview(_rightInput, cell, field, point, Kokkos::ALL(), Kokkos::ALL()) :
1356  Kokkos::subview(_rightInput, field, point, Kokkos::ALL(), Kokkos::ALL())):
1357  ( hasCell ? Kokkos::subview(_rightInput, cell, point, Kokkos::ALL(), Kokkos::ALL()) :
1358  Kokkos::subview(_rightInput, point, Kokkos::ALL(), Kokkos::ALL())));
1359 
1360  const ordinal_type iend = result.extent(0);
1361  const ordinal_type jend = result.extent(1);
1362 
1363  switch (leftRank) {
1364  case 4: {
1365  if (_isTranspose) {
1366  const size_type kend = right.extent(0);
1367  for (ordinal_type i=0; i<iend; ++i)
1368  for (ordinal_type j=0; j<jend; ++j) {
1369  result(i, j) = value_type(0);
1370  for (size_type k=0; k<kend; ++k)
1371  result(i, j) += left(k, i) * right(k, j);
1372  }
1373  } else {
1374  const size_type kend = right.extent(0);
1375  for (ordinal_type i=0; i<iend; ++i)
1376  for (ordinal_type j=0; j<jend; ++j) {
1377  result(i, j) = value_type(0);
1378  for (size_type k=0; k<kend; ++k)
1379  result(i, j) += left(i, k) * right(k, j);
1380  }
1381  }
1382  break;
1383  }
1384  case 3: { //matrix is diagonal
1385  for (ordinal_type i=0; i<iend; ++i)
1386  for (ordinal_type j=0; j<jend; ++j)
1387  result(i, j) = left(i) * right(i, j);
1388  break;
1389  }
1390  case 2: { //matrix is a scaled identity
1391  for (ordinal_type i=0; i<iend; ++i)
1392  for (ordinal_type j=0; j<jend; ++j)
1393  result(i, j) = left() * right(i, j);
1394  break;
1395  }
1396  }
1397  }
1398  };
1399  } //namespace
1400 
1401  template<typename DeviceType>
1402  template<typename outputValueType, class ...outputProperties,
1403  typename leftInputValueType, class ...leftInputProperties,
1404  typename rightInputValueType, class ...rightInputProperties>
1405  void
1406  ArrayTools<DeviceType>::Internal::
1407  matmatProduct( Kokkos::DynRankView<outputValueType, outputProperties...> output,
1408  const Kokkos::DynRankView<leftInputValueType, leftInputProperties...> leftInput,
1409  const Kokkos::DynRankView<rightInputValueType,rightInputProperties...> rightInput,
1410  const bool hasField,
1411  const bool isTranspose ) {
1412 
1413  typedef Kokkos::DynRankView<outputValueType, outputProperties...> OutputViewType;
1414  typedef const Kokkos::DynRankView<leftInputValueType, leftInputProperties...> leftInputViewType;
1415  typedef const Kokkos::DynRankView<rightInputValueType,rightInputProperties...> rightInputViewType;
1416  typedef FunctorArrayTools::F_matmatProduct<OutputViewType, leftInputViewType, rightInputViewType> FunctorType;
1417 
1418  const size_type loopSize = ( hasField ? output.extent(0)*output.extent(1)*output.extent(2) :
1419  output.extent(0)*output.extent(1) );
1420  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
1421  Kokkos::parallel_for( policy, FunctorType(output, leftInput, rightInput, hasField, isTranspose) );
1422  }
1423 
1424 
1425 
1426 
1427  template<typename DeviceType>
1428  template<typename outputFieldValueType, class ...outputFieldProperties,
1429  typename inputDataValueType, class ...inputDataProperties,
1430  typename inputFieldValueType, class ...inputFieldProperties>
1431  void
1432  ArrayTools<DeviceType>::
1433  matmatProductDataField( Kokkos::DynRankView<outputFieldValueType,outputFieldProperties...> outputFields,
1434  const Kokkos::DynRankView<inputDataValueType, inputDataProperties...> inputData,
1435  const Kokkos::DynRankView<inputFieldValueType, inputFieldProperties...> inputFields,
1436  const char transpose ) {
1437 
1438 #ifdef HAVE_INTREPID2_DEBUG
1439  {
1440  /*
1441  * Check array rank and spatial dimension range (if applicable)
1442  * (1) inputData(C,P), (C,P,D) or (C,P,D,D); P=1 is admissible to allow multiply by const. data
1443  * (2) inputFields(C,F,P,D,D) or (F,P,D,D);
1444  * (3) outputFields(C,F,P,D,D)
1445  */
1446  // (1) inputData is (C,P), (C, P, D) or (C, P, D, D) and 1 <= D <= 3 is required
1447  INTREPID2_TEST_FOR_EXCEPTION( inputData.rank() < 2 ||
1448  inputData.rank() > 4, std::invalid_argument,
1449  ">>> ERROR (ArrayTools::matmatProductDataField): inputData must have rank 2,3 or 4" );
1450  if (inputData.rank() > 2) {
1451  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(2) < 1 ||
1452  inputData.extent(2) > 3, std::invalid_argument,
1453  ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension (2) must be 1,2 or 3" );
1454  }
1455  if (inputData.rank() == 4) {
1456  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(3) < 1 ||
1457  inputData.extent(3) > 3, std::invalid_argument,
1458  ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension (3) must be 1,2 or 3" );
1459  }
1460 
1461  // (2) inputFields is (C,F,P,D,D) or (F,P,D,D) and 1 <= (dimension(rank-1), (rank-2)) <= 3 is required.
1462  INTREPID2_TEST_FOR_EXCEPTION( inputFields.rank() < 4 ||
1463  inputFields.rank() > 5, std::invalid_argument,
1464  ">>> ERROR (ArrayTools::matmatProductDataField): inputFields must have rank 4 or 5");
1465  INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(inputFields.rank()-1) < 1 ||
1466  inputFields.extent(inputFields.rank()-1) > 3, std::invalid_argument,
1467  ">>> ERROR (ArrayTools::matmatProductDataField): inputFields dimension (rank-1) must be 1,2 or 3" );
1468  INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(inputFields.rank()-2) < 1 ||
1469  inputFields.extent(inputFields.rank()-2) > 3, std::invalid_argument,
1470  ">>> ERROR (ArrayTools::matmatProductDataField): inputFields dimension (rank-2) must be 1,2 or 3" );
1471 
1472  // (3) outputFields is (C,F,P,D,D)
1473  INTREPID2_TEST_FOR_EXCEPTION( outputFields.rank() != 5, std::invalid_argument,
1474  ">>> ERROR (ArrayTools::matmatProductDataField): outputFields must have rank 5" );
1475  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(3) < 1 ||
1476  outputFields.extent(3) > 3, std::invalid_argument,
1477  ">>> ERROR (ArrayTools::matmatProductDataField): outputFields dimension (3) must be 1,2 or 3" );
1478  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(4) < 1 ||
1479  outputFields.extent(4) > 3, std::invalid_argument,
1480  ">>> ERROR (ArrayTools::matmatProductDataField): outputFields dimension (4) must be 1,2 or 3" );
1481 
1482  /*
1483  * Dimension cross-checks:
1484  * (1) inputData vs. inputFields
1485  * (2) outputFields vs. inputData
1486  * (3) outputFields vs. inputFields
1487  *
1488  * Cross-check (2): outputFields(C,F,P,D,D) vs. inputData(C,P), (C,P,D) or (C,P,D,D):
1489  * dimensions C, and D must match in all cases, dimension P must match only when non-constant
1490  * data is specified (P>1). Do not check P dimensions with constant data, i.e., when P=1 in
1491  * inputData(C,1,...)
1492  */
1493  if (inputData.extent(1) > 1) { // check P dimension if P>1 in inputData
1494  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(2) != inputData.extent(1), std::invalid_argument,
1495  ">>> ERROR (ArrayTools::matmatProductDataField): outputFields dimension (2) does not match to inputData dimension (1)" );
1496  }
1497  if (inputData.rank() == 2) { // inputData(C,P) -> C match
1498  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(0) != inputData.extent(0), std::invalid_argument,
1499  ">>> ERROR (ArrayTools::matmatProductDataField): outputFields dimension (0) does not match to inputData dimension (0)" );
1500  }
1501  if (inputData.rank() == 3) { // inputData(C,P,D) -> C, D match
1502  const size_type f1[] = { 0, 3, 4 }, f2[] = { 0, 2, 2 };
1503  for (size_type i=0; i<3; ++i) {
1504  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputData.extent(f2[i]), std::invalid_argument,
1505  ">>> ERROR (ArrayTools::matmatProductDataField): outputFields dimension does not match to inputData dimension" );
1506  }
1507  }
1508  if (inputData.rank() == 4) { // inputData(C,P,D,D) -> C, D, D match
1509  const size_type f1[] = { 0, 3, 4 }, f2[] = { 0, 2, 3 };
1510  for (size_type i=0; i<3; ++i) {
1511  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputData.extent(f2[i]), std::invalid_argument,
1512  ">>> ERROR (ArrayTools::matmatProductDataField): outputFields dimension does not match to inputData dimension" );
1513  }
1514  }
1515 
1516  /*
1517  * Cross-checks (1,3):
1518  */
1519  if (inputFields.rank() == 5) {
1520  // Cross-check (1): inputData(C,P), (C,P,D) or (C,P,D,D) vs. inputFields(C,F,P,D,D): dimensions C, P, D must match
1521  if (inputData.extent(1) > 1) { // check P dimension if P>1 in inputData
1522  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(1) != inputFields.extent(2), std::invalid_argument,
1523  ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension (1) does not match to inputFields dimension (2)" );
1524  }
1525  if (inputData.rank() == 2) { // inputData(C,P) -> C match
1526  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(0) != inputFields.extent(0), std::invalid_argument,
1527  ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension (0) does not match to inputFields dimension (0)" );
1528  }
1529  if (inputData.rank() == 3) { // inputData(C,P,D) -> C, D match
1530 
1531  const size_type f1[] = { 0, 2, 2 }, f2[] = { 0, 3, 4 };
1532  for (size_type i=0; i<3; ++i) {
1533  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
1534  ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension does not match to inputFields dimension" );
1535  }
1536  }
1537  if (inputData.rank() == 4) { // inputData(C,P,D,D) -> C, D, D match
1538  const size_type f1[] = { 0, 2, 3 }, f2[] = { 0, 3, 4 };
1539  for (size_type i=0; i<3; ++i) {
1540  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
1541  ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension does not match to inputFields dimension" );
1542  }
1543  }
1544 
1545  // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(C,F,P,D,D): all dimensions C, F, P, D must match
1546  for (size_type i=0; i<outputFields.rank(); ++i) {
1547  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(i) != inputFields.extent(i), std::invalid_argument,
1548  ">>> ERROR (ArrayTools::matmatProductDataField): outputFields dimension does not match to inputFields dimension" );
1549  }
1550  } else {
1551  // Cross-check (1): inputData(C,P), (C,P,D) or (C,P,D,D) vs. inputFields(F,P,D,D): dimensions P, D must match
1552  if (inputData.extent(1) > 1) { // check P if P>1 in inputData
1553  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(1) != inputFields.extent(1), std::invalid_argument,
1554  ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension (1) does not match to inputFields dimension (1)" );
1555  }
1556  if (inputData.rank() == 3) {
1557  const size_type f1[] = { 2, 2 }, f2[] = { 2, 3 };
1558  for (size_type i=0; i<2; ++i) {
1559  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
1560  ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension does not match to inputFields dimension" );
1561  }
1562  }
1563  if (inputData.rank() == 4) {
1564  const size_type f1[] = { 2, 3 }, f2[] = { 2, 3 };
1565  for (size_type i=0; i<2; ++i) {
1566  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
1567  ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension does not match to inputFields dimension" );
1568  }
1569  }
1570 
1571  // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(F,P,D,D): dimensions F, P, D must match
1572  {
1573  const size_type f1[] = { 1, 2, 3, 4 }, f2[] = { 0, 1, 2, 3 };
1574  for (size_type i=0; i<4; ++i) {
1575  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
1576  ">>> ERROR (ArrayTools::matmatProductDataField): outputFields dimension does not match to inputFields dimension" );
1577  }
1578  }
1579  }
1580  }
1581 #endif
1582  // body
1583  ArrayTools<DeviceType>::Internal::matmatProduct( outputFields,
1584  inputData,
1585  inputFields,
1586  true,
1587  transpose == 't' || transpose == 'T' );
1588  }
1589 
1590 
1591 
1592  template<typename DeviceType>
1593  template<typename outputDataValueType, class ...outputDataProperties,
1594  typename inputDataLeftValueType, class ...inputDataLeftProperties,
1595  typename inputDataRightValueType, class ...inputDataRightProperties>
1596  void
1597  ArrayTools<DeviceType>::matmatProductDataData( Kokkos::DynRankView<outputDataValueType, outputDataProperties...> outputData,
1598  const Kokkos::DynRankView<inputDataLeftValueType, inputDataLeftProperties...> inputDataLeft,
1599  const Kokkos::DynRankView<inputDataRightValueType,inputDataRightProperties...> inputDataRight,
1600  const char transpose ) {
1601 
1602 #ifdef HAVE_INTREPID2_DEBUG
1603  {
1604  /*
1605  * Check array rank and spatial dimension range (if applicable)
1606  * (1) inputDataLeft(C,P), (C,P,D) or (C,P,D,D); P=1 is admissible to allow multiply by const. left data
1607  * (2) inputDataRight(C,P,D,D) or (P,D,D);
1608  * (3) outputData(C,P,D,D)
1609  */
1610  // (1) inputDataLeft is (C,P), (C,P,D) or (C,P,D,D) and 1 <= D <= 3 is required
1611  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.rank() < 2 ||
1612  inputDataLeft.rank() > 4, std::invalid_argument,
1613  ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft must have rank 2,3 or 4" );
1614  if (inputDataLeft.rank() > 2) {
1615  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(2) < 1 ||
1616  inputDataLeft.extent(2) > 3, std::invalid_argument,
1617  ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft dimension (2) must be 1,2 or 3" );
1618  }
1619  if (inputDataLeft.rank() == 4) {
1620  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(3) < 1 ||
1621  inputDataLeft.extent(3) > 3, std::invalid_argument,
1622  ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft dimension (3) must be 1,2 or 3" );
1623  }
1624 
1625  // (2) inputDataRight is (C,P,D,D) or (P,D,D) and 1 <= (D=dimension(rank-1),(rank-2)) <= 3 is required.
1626  INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.rank() < 3 ||
1627  inputDataRight.rank() > 4, std::invalid_argument,
1628  ">>> ERROR (ArrayTools::matmatProductDataData): inputDataRight must have rank 3 or 4" );
1629  INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.extent(inputDataRight.rank()-1) < 1 ||
1630  inputDataRight.extent(inputDataRight.rank()-1) > 3, std::invalid_argument,
1631  ">>> ERROR (ArrayTools::matmatProductDataData): inputDataRight (rank-1) must be 1,2 or 3" );
1632  INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.extent(inputDataRight.rank()-2) < 1 ||
1633  inputDataRight.extent(inputDataRight.rank()-2) > 3, std::invalid_argument,
1634  ">>> ERROR (ArrayTools::matmatProductDataData): inputDataRight (rank-2) must be 1,2 or 3" );
1635 
1636  // (3) outputData is (C,P,D,D)
1637  INTREPID2_TEST_FOR_EXCEPTION( outputData.rank() != 4, std::invalid_argument,
1638  ">>> ERROR (ArrayTools::matmatProductDataData): outputData must have rank 4" );
1639  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(2) < 1 ||
1640  outputData.extent(2) > 3, std::invalid_argument,
1641  ">>> ERROR (ArrayTools::matmatProductDataData): outputData (2) must be 1,2 or 3" );
1642  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(3) < 1 ||
1643  outputData.extent(3) > 3, std::invalid_argument,
1644  ">>> ERROR (ArrayTools::matmatProductDataData): outputData (3) must be 1,2 or 3" );
1645 
1646  /*
1647  * Dimension cross-checks:
1648  * (1) inputDataLeft vs. inputDataRight
1649  * (2) outputData vs. inputDataLeft
1650  * (3) outputData vs. inputDataRight
1651  *
1652  * Cross-check (2): outputData(C,P,D,D) vs. inputDataLeft(C,P), (C,P,D) or (C,P,D,D):
1653  * dimensions C, and D must match in all cases, dimension P must match only when non-constant
1654  * data is specified (P>1). Do not check P dimensions with constant left data, i.e., when P=1 in
1655  * inputDataLeft(C,1,...)
1656  */
1657  if (inputDataLeft.extent(1) > 1) { // check P dimension if P>1 in inputDataLeft
1658  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(1) != inputDataLeft.extent(1), std::invalid_argument,
1659  ">>> ERROR (ArrayTools::matmatProductDataData): outputData dimension (1) does not match to inputDataLeft dimension (1)" );
1660  }
1661  if (inputDataLeft.rank() == 2) { // inputDataLeft(C,P): check C
1662  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(0) != inputDataLeft.extent(0), std::invalid_argument,
1663  ">>> ERROR (ArrayTools::matmatProductDataData): outputData dimension (0) does not match to inputDataLeft dimension (0)" );
1664  }
1665  if (inputDataLeft.rank() == 3) { // inputDataLeft(C,P,D): check C and D
1666  const size_type f1[] = { 0, 2, 3 }, f2[] = { 0, 2, 2 };
1667  for (size_type i=0; i<3; ++i) {
1668  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataLeft.extent(f2[i]), std::invalid_argument,
1669  ">>> ERROR (ArrayTools::matmatProductDataData): outputData dimension does not match to inputDataLeft dimension" );
1670  }
1671  }
1672  if (inputDataLeft.rank() == 4) { // inputDataLeft(C,P,D,D): check C and D
1673  const size_type f1[] = { 0, 2, 3 }, f2[] = { 0, 2, 3 };
1674  for (size_type i=0; i<3; ++i) {
1675  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataLeft.extent(f2[i]), std::invalid_argument,
1676  ">>> ERROR (ArrayTools::matmatProductDataData): outputData dimension does not match to inputDataLeft dimension" );
1677  }
1678  }
1679 
1680  /*
1681  * Cross-checks (1,3):
1682  */
1683  if (inputDataRight.rank() == 4) {
1684  // Cross-check (1): inputDataLeft(C,P), (C,P,D), or (C,P,D,D) vs. inputDataRight(C,P,D,D): dimensions C, P, D must match
1685  if (inputDataLeft.extent(1) > 1) { // check P dimension if P>1 in inputDataLeft
1686  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(1) != inputDataRight.extent(1), std::invalid_argument,
1687  ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft dimension (1) does not match to inputDataRight dimension (1)" );
1688  }
1689  if (inputDataLeft.rank() == 2) { // inputDataLeft(C,P): check C
1690  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(0) != inputDataRight.extent(0), std::invalid_argument,
1691  ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft dimension (0) does not match to inputDataRight dimension (0)" );
1692  }
1693  if (inputDataLeft.rank() == 3) { // inputDataLeft(C,P,D): check C and D
1694  const size_type f1[] = { 0, 2, 2 }, f2[] = { 0, 2, 3 };
1695  for (size_type i=0; i<3; ++i) {
1696  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1697  ">>> ERROR (ArrayTools::matmatProductDataData): intputDataLeft dimension does not match to inputDataRight dimension" );
1698  }
1699  }
1700  if (inputDataLeft.rank() == 4) { // inputDataLeft(C,P,D,D): check C and D
1701  const size_type f1[] = { 0, 2, 3 }, f2[] = { 0, 2, 3 };
1702  for (size_type i=0; i<3; ++i) {
1703  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1704  ">>> ERROR (ArrayTools::matmatProductDataData): intputDataLeft dimension does not match to inputDataRight dimension" );
1705  }
1706  }
1707 
1708  // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(C,P,D,D): all dimensions C, P, D must match
1709  for (size_type i=0; i< outputData.rank(); ++i) {
1710  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(i) != inputDataRight.extent(i), std::invalid_argument,
1711  ">>> ERROR (ArrayTools::matmatProductDataData): outputData dimension does not match to inputDataRight dimension" );
1712  }
1713  } else {
1714  // Cross-check (1): inputDataLeft(C,P), (C,P,D), or (C,P,D,D) vs. inputDataRight(P,D,D): dimensions P, D must match
1715  if (inputDataLeft.extent(1) > 1) { // check P if P>1 in inputData
1716  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(1) != inputDataRight.extent(0), std::invalid_argument,
1717  ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft dimension (1) does not match to inputDataRight dimension (0)" );
1718  }
1719  if (inputDataLeft.rank() == 3) { // inputDataLeft(C,P,D): check D
1720  const size_type f1[] = { 2, 2 }, f2[] = { 1, 2 };
1721  for (size_type i=0; i<2; ++i) {
1722  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1723  ">>> ERROR (ArrayTools::matmatProductDataData): intputDataLeft dimension does not match to inputDataRight dimension" );
1724  }
1725  }
1726  if (inputDataLeft.rank() == 4) { // inputDataLeft(C,P,D,D): check D
1727  const size_type f1[] = { 2, 3 }, f2[] = { 1, 2 };
1728  for (size_type i=0; i<2; ++i) {
1729  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1730  ">>> ERROR (ArrayTools::matmatProductDataData): intputDataLeft dimension does not match to inputDataRight dimension" );
1731  }
1732  }
1733  // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(P,D,D): dimensions P, D must match
1734  {
1735  const size_type f1[] = { 1, 2, 3 }, f2[] = { 0, 1, 2 };
1736  for (size_type i=0; i<3; ++i) {
1737  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1738  ">>> ERROR (ArrayTools::matmatProductDataData): outputData dimension does not match to inputDataRight dimension" );
1739  }
1740  }
1741  }
1742  }
1743 #endif
1744  // body
1745  ArrayTools<DeviceType>::Internal::matmatProduct( outputData,
1746  inputDataLeft,
1747  inputDataRight,
1748  false,
1749  transpose == 't' || transpose == 'T' );
1750  }
1751 
1752 } // end namespace Intrepid2
1753 #endif
Intrepid2::ArrayTools::matmatProductDataData
static void matmatProductDataData(Kokkos::DynRankView< outputDataValueType, outputDataProperties...> outputData, const Kokkos::DynRankView< inputDataLeftValueType, inputDataLeftProperties...> inputDataLeft, const Kokkos::DynRankView< inputDataRightValueType, inputDataRightProperties...> inputDataRight, const char transpose= 'N')
There are two use cases: (1) matrix-matrix product of a rank-4 container inputDataRight with dimensio...
Definition: Intrepid2_ArrayToolsDefTensor.hpp:1597
Intrepid2::FunctorArrayTools::F_crossProduct
Functor for crossProduct see Intrepid2::ArrayTools for more.
Definition: Intrepid2_ArrayToolsDefTensor.hpp:26
Intrepid2::FunctorArrayTools::F_outerProduct
Functor for outerProduct see Intrepid2::ArrayTools for more.
Definition: Intrepid2_ArrayToolsDefTensor.hpp:352
Intrepid2::ArrayTools::crossProductDataData
static void crossProductDataData(Kokkos::DynRankView< outputDataValueType, outputDataProperties...> outputData, const Kokkos::DynRankView< inputDataLeftValueType, inputDataLeftProperties...> inputDataLeft, const Kokkos::DynRankView< inputDataRightValueType, inputDataRightProperties...> inputDataRight)
There are two use cases: (1) cross product of a rank-3 container inputDataRight with dimensions (C...
Definition: Intrepid2_ArrayToolsDefTensor.hpp:235
Intrepid2::ArrayTools::crossProductDataField
static void crossProductDataField(Kokkos::DynRankView< outputFieldValueType, outputFieldProperties...> outputFields, const Kokkos::DynRankView< inputDataValueType, inputDataProperties...> inputData, const Kokkos::DynRankView< inputFieldValueType, inputFieldProperties...> inputFields)
There are two use cases: (1) cross product of a rank-4 container inputFields with dimensions (C...
Definition: Intrepid2_ArrayToolsDefTensor.hpp:112
Intrepid2::FunctorArrayTools::F_matmatProduct
Functor for matmatProduct see Intrepid2::ArrayTools for more.
Definition: Intrepid2_ArrayToolsDefTensor.hpp:1311
Intrepid2::ArrayTools::matvecProductDataField
static void matvecProductDataField(Kokkos::DynRankView< outputFieldValueType, outputFieldProperties...> outputFields, const Kokkos::DynRankView< inputDataValueType, inputDataProperties...> inputData, const Kokkos::DynRankView< inputFieldValueType, inputFieldProperties...> inputFields, const char transpose= 'N')
There are two use cases: (1) matrix-vector product of a rank-4 container inputFields with dimensions ...
Definition: Intrepid2_ArrayToolsDefTensor.hpp:1001
Intrepid2::FunctorArrayTools::F_matvecProduct
Functor for matvecProduct; new version avoids both subviews and branching. See Intrepid2::ArrayTools ...
Definition: Intrepid2_ArrayToolsDefTensor.hpp:648
Intrepid2::ArrayTools::outerProductDataField
static void outerProductDataField(Kokkos::DynRankView< outputFieldValueType, outputFieldProperties...> outputFields, const Kokkos::DynRankView< inputDataValueType, inputDataProperties...> inputData, const Kokkos::DynRankView< inputFieldValueType, inputFieldProperties...> inputFields)
There are two use cases: (1) outer product of a rank-4 container inputFields with dimensions (C...
Definition: Intrepid2_ArrayToolsDefTensor.hpp:432
Intrepid2::ArrayTools::matmatProductDataField
static void matmatProductDataField(Kokkos::DynRankView< outputFieldValueType, outputFieldProperties...> outputFields, const Kokkos::DynRankView< inputDataValueType, inputDataProperties...> inputData, const Kokkos::DynRankView< inputFieldValueType, inputFieldProperties...> inputFields, const char transpose= 'N')
There are two use cases: (1) matrix-matrix product of a rank-5 container inputFields with dimensions ...
Definition: Intrepid2_ArrayToolsDefTensor.hpp:1433
Intrepid2::ArrayTools::outerProductDataData
static void outerProductDataData(Kokkos::DynRankView< outputDataValueType, outputDataProperties...> outputData, const Kokkos::DynRankView< inputDataLeftValuetype, inputDataLeftProperties...> inputDataLeft, const Kokkos::DynRankView< inputDataRightValueType, inputDataRightProperties...> inputDataRight)
There are two use cases: (1) outer product of a rank-3 container inputDataRight with dimensions (C...
Definition: Intrepid2_ArrayToolsDefTensor.hpp:539
Intrepid2::ArrayTools
Utility class that provides methods for higher-order algebraic manipulation of user-defined arrays...
Definition: Intrepid2_ArrayTools.hpp:44
Intrepid2::ArrayTools::matvecProductDataData
static void matvecProductDataData(Kokkos::DynRankView< outputDataValueType, outputDataProperties...> outputData, const Kokkos::DynRankView< inputDataLeftValueType, inputDataLeftProperties...> inputDataLeft, const Kokkos::DynRankView< inputDataRightValueType, inputDataRightProperties...> inputDataRight, const char transpose= 'N')
There are two use cases: (1) matrix-vector product of a rank-3 container inputDataRight with dimensio...
Definition: Intrepid2_ArrayToolsDefTensor.hpp:1155
Generated by   doxygen 1.8.5