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