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