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 >
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
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;
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
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
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
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
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
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;
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
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
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
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
664  inputDataLeft,
665  inputDataRight,
666  false );
667  }
668 
669 
670  namespace FunctorArrayTools {
674  template < typename OutputViewType,
675  typename leftInputViewType,
676  typename rightInputViewType>
678  OutputViewType _output;
679  const leftInputViewType _leftInput;
680  const rightInputViewType _rightInput;
681 
682  const bool _isTranspose;
683 
684  KOKKOS_INLINE_FUNCTION
685  F_matvecProduct(OutputViewType output_,
686  leftInputViewType leftInput_,
687  rightInputViewType rightInput_,
688  const bool isTranspose_)
689  : _output(output_), _leftInput(leftInput_), _rightInput(rightInput_), _isTranspose(isTranspose_) {}
690 
691  template<typename resultViewType,
692  typename leftViewType,
693  typename rightViewType>
694  KOKKOS_FORCEINLINE_FUNCTION
695  static void
696  apply_matvec_product( resultViewType &result,
697  const leftViewType &left,
698  const rightViewType &right,
699  const bool isTranspose) {
700  const ordinal_type iend = result.extent(0);
701  const ordinal_type jend = right.extent(0);
702 
703  typedef typename resultViewType::value_type value_type;
704 
705  switch (left.rank()) {
706  case 2:
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 += left(j, i)*right(j);
712  result(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 += left(i, j)*right(j);
719  result(i) = tmp;
720  }
721  }
722  break;
723  case 1: { //matrix is diagonal
724  for (ordinal_type i=0;i<iend;++i)
725  result(i) = left(i)*right(i);
726  break;
727  }
728  case 0: { //matrix is a scaled identity
729  const value_type val = left();
730  for (ordinal_type i=0;i<iend;++i) {
731  result(i) = val*right(i);
732  }
733  break;
734  }
735  }
736  }
737 
738  KOKKOS_INLINE_FUNCTION
739  void operator()(const ordinal_type cl,
740  const ordinal_type pt) const {
741  const auto rightRank = _rightInput.rank();
742  const auto leftRank = _leftInput.rank();
743 
744  auto result = Kokkos::subview(_output, cl, pt, Kokkos::ALL());
745 
746  const auto lpt = (_leftInput.extent(1) == 1 ? size_type(0) : pt);
747  const auto left = ( leftRank == 4 ? Kokkos::subview(_leftInput, cl, lpt, Kokkos::ALL(), Kokkos::ALL()) :
748  leftRank == 3 ? Kokkos::subview(_leftInput, cl, lpt, Kokkos::ALL()) :
749  Kokkos::subview(_leftInput, cl, lpt));
750 
751  const auto right = ( rightRank == 2 ? Kokkos::subview(_rightInput, pt, Kokkos::ALL()) :
752  Kokkos::subview(_rightInput, cl, pt, Kokkos::ALL()) );
753  apply_matvec_product( result, left, right, _isTranspose );
754  }
755 
756  KOKKOS_INLINE_FUNCTION
757  void operator()(const ordinal_type cl,
758  const ordinal_type bf,
759  const ordinal_type pt) const {
760  const auto rightRank = _rightInput.rank();
761  const auto leftRank = _leftInput.rank();
762 
763  auto result = Kokkos::subview(_output, cl, bf, pt, Kokkos::ALL());
764 
765  const auto lpt = (_leftInput.extent(1) == 1 ? size_type(0) : pt);
766  const auto left = ( leftRank == 4 ? Kokkos::subview(_leftInput, cl, lpt, Kokkos::ALL(), Kokkos::ALL()) :
767  leftRank == 3 ? Kokkos::subview(_leftInput, cl, lpt, Kokkos::ALL()) :
768  Kokkos::subview(_leftInput, cl, lpt));
769 
770  const auto right = ( rightRank == 3 ? Kokkos::subview(_rightInput, bf, pt, Kokkos::ALL()) :
771  Kokkos::subview(_rightInput, cl, bf, pt, Kokkos::ALL()));
772 
773  apply_matvec_product( result, left, right, _isTranspose );
774  }
775  };
776  } //namespace
777 
778  template<typename DeviceType>
779  template<typename outputValueType, class ...outputProperties,
780  typename leftInputValueType, class ...leftInputProperties,
781  typename rightInputValueType, class ...rightInputProperties>
782  void
784  matvecProduct( Kokkos::DynRankView<outputValueType, outputProperties...> output,
785  const Kokkos::DynRankView<leftInputValueType, leftInputProperties...> leftInput,
786  const Kokkos::DynRankView<rightInputValueType,rightInputProperties...> rightInput,
787  const bool hasField,
788  const bool isTranspose ) {
789 
790  typedef Kokkos::DynRankView<outputValueType, outputProperties...> OutputViewType;
791  typedef const Kokkos::DynRankView<leftInputValueType, leftInputProperties...> leftInputViewType;
792  typedef const Kokkos::DynRankView<rightInputValueType,rightInputProperties...> rightInputViewType;
794 
795  if (hasField) {
796  using range_policy_type = Kokkos::Experimental::MDRangePolicy
797  < ExecSpaceType, Kokkos::Experimental::Rank<3>, Kokkos::IndexType<ordinal_type> >;
798  range_policy_type policy( { 0, 0, 0 },
799  { output.extent(0), output.extent(1), output.extent(2) } );
800  Kokkos::parallel_for( policy, FunctorType(output, leftInput, rightInput, isTranspose) );
801  } else {
802  using range_policy_type = Kokkos::Experimental::MDRangePolicy
803  < ExecSpaceType, Kokkos::Experimental::Rank<2>, Kokkos::IndexType<ordinal_type> >;
804  range_policy_type policy( { 0, 0 },
805  { output.extent(0), output.extent(1) } );
806  Kokkos::parallel_for( policy, FunctorType(output, leftInput, rightInput, isTranspose) );
807  }
808  }
809 
810  template<typename DeviceType>
811  template<typename outputFieldValueType, class ...outputFieldProperties,
812  typename inputDataValueType, class ...inputDataProperties,
813  typename inputFieldValueType, class ...inputFieldProperties>
814  void
815  ArrayTools<DeviceType>::matvecProductDataField( Kokkos::DynRankView<outputFieldValueType,outputFieldProperties...> outputFields,
816  const Kokkos::DynRankView<inputDataValueType, inputDataProperties...> inputData,
817  const Kokkos::DynRankView<inputFieldValueType, inputFieldProperties...> inputFields,
818  const char transpose ) {
819 
820 #ifdef HAVE_INTREPID2_DEBUG
821  {
822  /*
823  * Check array rank and spatial dimension range (if applicable)
824  * (1) inputData(C,P), (C,P,D) or (C,P,D,D); P=1 is admissible to allow multiply by const. data
825  * (2) inputFields(C,F,P,D) or (F,P,D);
826  * (3) outputFields(C,F,P,D)
827  */
828  // (1) inputData is (C,P), (C, P, D) or (C, P, D, D) and 1 <= D <= 3 is required
829  INTREPID2_TEST_FOR_EXCEPTION( inputData.rank() < 2 || inputData.rank() > 4, std::invalid_argument,
830  ">>> ERROR (ArrayTools::matvecProductDataField): inputData must have rank 2,3 or 4" );
831  if (inputData.rank() > 2) {
832  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(2) < 1 ||
833  inputData.extent(2) > 3, std::invalid_argument,
834  ">>> ERROR (ArrayTools::matvecProductDataField): inputData dimension(2) must be 1,2 or 3");
835  }
836  if (inputData.rank() == 4) {
837  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(3) < 1 ||
838  inputData.extent(3) > 3, std::invalid_argument,
839  ">>> ERROR (ArrayTools::matvecProductDataField): inputData dimension(3) must be 1,2 or 3");
840  }
841 
842  // (2) inputFields is (C, F, P, D) or (F, P, D) and 1 <= (D=dimension(rank - 1)) <= 3 is required.
843  INTREPID2_TEST_FOR_EXCEPTION( inputFields.rank() < 3 ||
844  inputFields.rank() > 4, std::invalid_argument,
845  ">>> ERROR (ArrayTools::matvecProductDataField): inputFields must have rank 3 or 4" );
846  INTREPID2_TEST_FOR_EXCEPTION( (inputFields.rank()-1) < 1 ||
847  (inputFields.rank()-1) > 3, std::invalid_argument,
848  ">>> ERROR (ArrayTools::matvecProductDataField): inputFields dimension(rank-1) must be 1,2, or 3" );
849 
850  // (3) outputFields is (C,F,P,D)
851  INTREPID2_TEST_FOR_EXCEPTION( outputFields.rank() != 4, std::invalid_argument,
852  ">>> ERROR (ArrayTools::matvecProductDataField): outputFields must have rank 4" );
853  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(3) < 1 ||
854  outputFields.extent(3) > 3, std::invalid_argument,
855  ">>> ERROR (ArrayTools::matvecProductDataField): outputFields dimension(3) must be 1,2 or 3" );
856 
857  /*
858  * Dimension cross-checks:
859  * (1) inputData vs. inputFields
860  * (2) outputFields vs. inputData
861  * (3) outputFields vs. inputFields
862  *
863  * Cross-check (2): outputFields(C,F,P,D) vs. inputData(C,P), (C,P,D) or (C,P,D,D):
864  * dimensions C, and D must match in all cases, dimension P must match only when non-constant
865  * data is specified (P>1). Do not check P dimensions with constant data, i.e., when P=1 in
866  * inputData(C,1,...)
867  */
868  if (inputData.extent(1) > 1) { // check P dimension if P>1 in inputData
869  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(2) != inputData.extent(1), std::invalid_argument,
870  ">>> ERROR (ArrayTools::matvecProductDataField): outputFields dimension(2) must match to inputData dimension(1)" );
871  }
872  if (inputData.rank() == 2) { // inputData(C,P) -> C match
873  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(0) != inputData.extent(0), std::invalid_argument,
874  ">>> ERROR (ArrayTools::matvecProductDataField): outputFields dimension(0) must match to inputData dimension(0)" );
875  }
876  if (inputData.rank() == 3) { // inputData(C,P,D) -> C, D match
877  const size_type f1[] = { 0, 3 }, f2[] = { 0, 2 };
878  for (size_type i=0; i<2; ++i) {
879  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputData.extent(f2[i]), std::invalid_argument,
880  ">>> ERROR (ArrayTools::matvecProductDataField): outputFields dimension does not match to inputData dimension" );
881  }
882  }
883  if (inputData.rank() == 4) { // inputData(C,P,D,D) -> C, D, D match
884  const size_type f1[] = { 0, 3, 3 }, f2[] = { 0, 2, 3 };
885  for (size_type i=0; i<3; ++i) {
886  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputData.extent(f2[i]), std::invalid_argument,
887  ">>> ERROR (ArrayTools::matvecProductDataField): outputFields dimension does not match to inputData dimension" );
888  }
889  }
890 
891  /*
892  * Cross-checks (1,3):
893  */
894  if (inputFields.rank() == 4) {
895  // 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
896  if (inputData.extent(1) > 1) { // check P dimension if P>1 in inputData
897  INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(2) != inputData.extent(1), std::invalid_argument,
898  ">>> ERROR (ArrayTools::matvecProductDataField): inputFields dimension (2) does not match to inputData dimension (1)" );
899  }
900  if (inputData.rank() == 2) { // inputData(C,P) -> C match
901  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(0) != inputFields.extent(0), std::invalid_argument,
902  ">>> ERROR (ArrayTools::matvecProductDataField): inputData dimension (0) does not match to inputFields dimension (0)" );
903  }
904  if (inputData.rank() == 3) { // inputData(C,P,D) -> C, D match
905  const size_type f1[] = { 0, 2 }, f2[] = { 0, 3 };
906  for (size_type i=0; i<2; ++i) {
907  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
908  ">>> ERROR (ArrayTools::matvecProductDataField): inputData dimension does not match to inputFields dimension" );
909  }
910  }
911  if (inputData.rank() == 4) { // inputData(C,P,D,D) -> C, D, D match
912  const size_type f1[] = { 0, 2, 3 }, f2[] = { 0, 3, 3 };
913  for (size_type i=0; i<3; ++i) {
914  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
915  ">>> ERROR (ArrayTools::matvecProductDataField): inputData dimension does not match to inputFields dimension" );
916  }
917  }
918 
919  // Cross-check (3): outputFields(C,F,P,D) vs. inputFields(C,F,P,D): all dimensions C, F, P, D must match
920  for (size_type i=0; i<outputFields.rank(); ++i) {
921  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(i) != inputFields.extent(i), std::invalid_argument,
922  ">>> ERROR (ArrayTools::matvecProductDataField): outputFields dimension does not match to inputFields dimension" );
923  }
924  } else {
925  // Cross-check (1): inputData(C,P), (C,P,D) or (C,P,D,D) vs. inputFields(F,P,D): dimensions P, D must match
926  if (inputData.extent(1) > 1) { // check P if P>1 in inputData
927  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(1) != inputFields.extent(1), std::invalid_argument,
928  ">>> ERROR (ArrayTools::matvecProductDataField): inputData dimension(1) does not match to inputFields dimension (1)" );
929  }
930  if (inputData.rank() == 3) {
931  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(2) != inputFields.extent(2), std::invalid_argument,
932  ">>> ERROR (ArrayTools::matvecProductDataField): inputData dimension(2) does not match to inputFields dimension (2)" );
933  }
934  if (inputData.rank() == 4) {
935  const size_type f1[] = { 2, 3 }, f2[] = { 2, 2 };
936  for (size_type i=0; i<2; ++i) {
937  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
938  ">>> ERROR (ArrayTools::matvecProductDataField): inputData dimension does not match to inputFields dimension" );
939  }
940  }
941 
942  // Cross-check (3): outputFields(C,F,P,D) vs. inputFields(F,P,D): dimensions F, P, D must match
943  {
944  const size_type f1[] = { 1, 2, 3 }, f2[] = { 0, 1, 2 };
945  for (size_type i=0; i<3; ++i) {
946  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
947  ">>> ERROR (ArrayTools::matvecProductDataField): outputFields dimension does not match to inputFields dimension" );
948  }
949  }
950  }
951  }
952 #endif
953  // body
955  inputData,
956  inputFields,
957  true,
958  transpose == 't' || transpose == 'T' );
959  }
960 
961 
962 
963  template<typename DeviceType>
964  template<typename outputDataValueType, class ...outputDataProperties,
965  typename inputDataLeftValueType, class ...inputDataLeftProperties,
966  typename inputDataRightValueType, class ...inputDataRightProperties>
967  void
969  matvecProductDataData( Kokkos::DynRankView<outputDataValueType, outputDataProperties...> outputData,
970  const Kokkos::DynRankView<inputDataLeftValueType, inputDataLeftProperties...> inputDataLeft,
971  const Kokkos::DynRankView<inputDataRightValueType,inputDataRightProperties...> inputDataRight,
972  const char transpose ) {
973 
974 #ifdef HAVE_INTREPID2_DEBUG
975  {
976  /*
977  * Check array rank and spatial dimension range (if applicable)
978  * (1) inputDataLeft(C,P), (C,P,D) or (C,P,D,D); P=1 is admissible to allow multiply by const. left data
979  * (2) inputDataRight(C,P,D) or (P,D);
980  * (3) outputData(C,P,D)
981  */
982  // (1) inputDataLeft is (C,P), (C,P,D) or (C,P,D,D) and 1 <= D <= 3 is required
983  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.rank() < 2 ||
984  inputDataLeft.rank() > 4, std::invalid_argument,
985  ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft must have rank 2,3 or 4" );
986 
987  if (inputDataLeft.rank() > 2) {
988  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(2) < 1 ||
989  inputDataLeft.extent(2) > 3, std::invalid_argument,
990  ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension(2) must be 1, 2 or 3");
991  }
992  if (inputDataLeft.rank() == 4) {
993  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(3) < 1 ||
994  inputDataLeft.extent(3) > 3, std::invalid_argument,
995  ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension(3) must be 1, 2 or 3");
996  }
997 
998  // (2) inputDataRight is (C, P, D) or (P, D) and 1 <= (D=dimension(rank - 1)) <= 3 is required.
999  INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.rank() < 2 ||
1000  inputDataRight.rank() > 3, std::invalid_argument,
1001  ">>> ERROR (ArrayTools::matvecProductDataData): inputDataRight must have rank 2 or 3" );
1002  INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.extent(inputDataRight.rank()-1) < 1 ||
1003  inputDataRight.extent(inputDataRight.rank()-1) > 3, std::invalid_argument,
1004  ">>> ERROR (ArrayTools::matvecProductDataData): inputDataRight dimension (rank-1) must be 1,2 or 3" );
1005 
1006  // (3) outputData is (C,P,D)
1007  INTREPID2_TEST_FOR_EXCEPTION( outputData.rank() != 3, std::invalid_argument,
1008  ">>> ERROR (ArrayTools::matvecProductDataData): outputData must have rank 3" );
1009  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(2) < 1 ||
1010  outputData.extent(2) > 3, std::invalid_argument,
1011  ">>> ERROR (ArrayTools::matvecProductDataData): outputData dimension(2) must be 1, 2 or 3");
1012 
1013  /*
1014  * Dimension cross-checks:
1015  * (1) inputDataLeft vs. inputDataRight
1016  * (2) outputData vs. inputDataLeft
1017  * (3) outputData vs. inputDataRight
1018  *
1019  * Cross-check (2): outputData(C,P,D) vs. inputDataLeft(C,P), (C,P,D) or (C,P,D,D):
1020  * dimensions C, and D must match in all cases, dimension P must match only when non-constant
1021  * data is specified (P>1). Do not check P dimensions with constant left data, i.e., when P=1 in
1022  * inputDataLeft(C,1,...)
1023  */
1024  if (inputDataLeft.extent(1) > 1) { // check P dimension if P>1 in inputDataLeft
1025  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(1) != inputDataLeft.extent(1), std::invalid_argument,
1026  ">>> ERROR (ArrayTools::matvecProductDataData): outputData dimension(1) muat match inputDataLeft dimension(1)" );
1027  }
1028  if (inputDataLeft.rank() == 2) { // inputDataLeft(C,P): check C
1029  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(0) != inputDataLeft.extent(0), std::invalid_argument,
1030  ">>> ERROR (ArrayTools::matvecProductDataData): outputData dimension(0) muat match inputDataLeft dimension(0)" );
1031  }
1032  if (inputDataLeft.rank() == 3) { // inputDataLeft(C,P,D): check C and D
1033  const size_type f1[] = { 0, 2 }, f2[] = { 0, 2 };
1034  for (size_type i=0; i<2; ++i) {
1035  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataLeft.extent(f2[i]), std::invalid_argument,
1036  ">>> ERROR (ArrayTools::matvecProductDataData): outputData dimension muat match inputDataLeft dimension" );
1037  }
1038  }
1039  if (inputDataLeft.rank() == 4) { // inputDataLeft(C,P,D,D): check C and D
1040  const size_type f1[] = { 0, 2, 2 }, f2[] = { 0, 2, 3 };
1041  for (size_type i=0; i<3; ++i) {
1042  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataLeft.extent(f2[i]), std::invalid_argument,
1043  ">>> ERROR (ArrayTools::matvecProductDataData): outputData dimension muat match inputDataLeft dimension" );
1044  }
1045  }
1046 
1047  /*
1048  * Cross-checks (1,3):
1049  */
1050  if (inputDataRight.rank() == 3) {
1051  // 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
1052  if (inputDataLeft.extent(1) > 1) { // check P dimension if P>1 in inputDataLeft
1053  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(1) != inputDataRight.extent(1), std::invalid_argument,
1054  ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension (1) muat match to inputDataRight dimension (1)" );
1055  }
1056  if (inputDataLeft.rank() == 2) { // inputDataLeft(C,P): check C
1057  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(0) != inputDataRight.extent(0), std::invalid_argument,
1058  ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension (0) muat match to inputDataRight dimension (0)" );
1059  }
1060  if (inputDataLeft.rank() == 3) { // inputDataLeft(C,P,D): check C and D
1061  const size_type f1[] = { 0, 2 }, f2[] = { 0, 2 };
1062  for (size_type i=0; i<2; ++i) {
1063  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1064  ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension muat match to inputDataRight dimension" );
1065  }
1066  }
1067  if (inputDataLeft.rank() == 4) { // inputDataLeft(C,P,D,D): check C and D
1068  const size_type f1[] = { 0, 2, 3 }, f2[] = { 0, 2, 2 };
1069  for (size_type i=0; i<3; ++i) {
1070  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1071  ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension muat match to inputDataRight dimension" );
1072  }
1073  }
1074 
1075  // Cross-check (3): outputData(C,P,D) vs. inputDataRight(C,P,D): all dimensions C, P, D must match
1076  for (size_type i=0; i<outputData.rank(); ++i) {
1077  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(i) != inputDataRight.extent(i), std::invalid_argument,
1078  ">>> ERROR (ArrayTools::matvecProductDataData): outputData dimension muat match to inputDataRight dimension" );
1079  }
1080  } else {
1081  // Cross-check (1): inputDataLeft(C,P), (C,P,D), or (C,P,D,D) vs. inputDataRight(P,D): dimensions P, D must match
1082  if (inputDataLeft.extent(1) > 1) { // check P if P>1 in inputData
1083  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(1) != inputDataRight.extent(0), std::invalid_argument,
1084  ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension (1) does mot match to inputDataright dimension (0)" );
1085  }
1086  if (inputDataLeft.rank() == 3) { // inputDataLeft(C,P,D): check D
1087  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(2) != inputDataRight.extent(1), std::invalid_argument,
1088  ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension (2) does mot match to inputDataright dimension (1)" );
1089  }
1090  if (inputDataLeft.rank() == 4) { // inputDataLeft(C,P,D,D): check D
1091  const size_type f1[] = { 2, 3 }, f2[] = { 1, 1 };
1092  for (size_type i=0; i<2; ++i) {
1093  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1094  ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension muat match to inputDataRight dimension" );
1095  }
1096  }
1097 
1098  // Cross-check (3): outputData(C,P,D) vs. inputDataRight(P,D): dimensions P, D must match
1099  {
1100  const size_type f1[] = { 1, 2 }, f2[] = { 0, 1 };
1101  for (size_type i=0; i<2; ++i) {
1102  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1103  ">>> ERROR (ArrayTools::matvecProductDataData): outputData dimension muat match to inputDataRight dimension" );
1104  }
1105  }
1106  }
1107  }
1108 #endif
1109  // body
1111  inputDataLeft,
1112  inputDataRight,
1113  false,
1114  transpose == 't' || transpose == 'T' );
1115  }
1116 
1117 
1118  namespace FunctorArrayTools {
1122  template < typename OutputViewType, typename leftInputViewType, typename rightInputViewType >
1124  OutputViewType _output;
1125  leftInputViewType _leftInput;
1126  rightInputViewType _rightInput;
1127  const bool _hasField, _isTranspose;
1128  typedef typename leftInputViewType::value_type value_type;
1129 
1130  KOKKOS_INLINE_FUNCTION
1131  F_matmatProduct(OutputViewType output_,
1132  leftInputViewType leftInput_,
1133  rightInputViewType rightInput_,
1134  const bool hasField_,
1135  const bool isTranspose_)
1136  : _output(output_), _leftInput(leftInput_), _rightInput(rightInput_),
1137  _hasField(hasField_), _isTranspose(isTranspose_) {}
1138 
1139  KOKKOS_INLINE_FUNCTION
1140  void operator()(const size_type iter) const {
1141  size_type cell(0), field(0), point(0);
1142  size_type leftRank = _leftInput.rank();
1143  size_type rightRank = _rightInput.rank();
1144 
1145  if (_hasField)
1146  unrollIndex( cell, field, point,
1147  _output.extent(0),
1148  _output.extent(1),
1149  _output.extent(2),
1150  iter );
1151  else
1152  unrollIndex( cell, point,
1153  _output.extent(0),
1154  _output.extent(1),
1155  iter );
1156 
1157  auto result = ( _hasField ? Kokkos::subview(_output, cell, field, point, Kokkos::ALL(), Kokkos::ALL()) :
1158  Kokkos::subview(_output, cell, point, Kokkos::ALL(), Kokkos::ALL()));
1159 
1160  const auto lpoint = (_leftInput.extent(1) == 1 ? 0 : point);
1161  auto left = ( leftRank == 4 ? Kokkos::subview(_leftInput, cell, lpoint, Kokkos::ALL(), Kokkos::ALL()) :
1162  leftRank == 3 ? Kokkos::subview(_leftInput, cell, lpoint, Kokkos::ALL()) :
1163  Kokkos::subview(_leftInput, cell, lpoint) );
1164 
1165  //temporary
1166  const bool hasCell = (_hasField ? rightRank == 5 : rightRank == 4);
1167  auto right = ( _hasField ? ( hasCell ? Kokkos::subview(_rightInput, cell, field, point, Kokkos::ALL(), Kokkos::ALL()) :
1168  Kokkos::subview(_rightInput, field, point, Kokkos::ALL(), Kokkos::ALL())):
1169  ( hasCell ? Kokkos::subview(_rightInput, cell, point, Kokkos::ALL(), Kokkos::ALL()) :
1170  Kokkos::subview(_rightInput, point, Kokkos::ALL(), Kokkos::ALL())));
1171 
1172  const ordinal_type iend = result.extent(0);
1173  const ordinal_type jend = result.extent(1);
1174 
1175  switch (leftRank) {
1176  case 4: {
1177  if (_isTranspose) {
1178  const size_type kend = right.extent(0);
1179  for (ordinal_type i=0; i<iend; ++i)
1180  for (ordinal_type j=0; j<jend; ++j) {
1181  result(i, j) = value_type(0);
1182  for (size_type k=0; k<kend; ++k)
1183  result(i, j) += left(k, i) * right(k, j);
1184  }
1185  } else {
1186  const size_type kend = right.extent(0);
1187  for (ordinal_type i=0; i<iend; ++i)
1188  for (ordinal_type j=0; j<jend; ++j) {
1189  result(i, j) = value_type(0);
1190  for (size_type k=0; k<kend; ++k)
1191  result(i, j) += left(i, k) * right(k, j);
1192  }
1193  }
1194  break;
1195  }
1196  case 3: { //matrix is diagonal
1197  for (ordinal_type i=0; i<iend; ++i)
1198  for (ordinal_type j=0; j<jend; ++j)
1199  result(i, j) = left(i) * right(i, j);
1200  break;
1201  }
1202  case 2: { //matrix is a scaled identity
1203  for (ordinal_type i=0; i<iend; ++i)
1204  for (ordinal_type j=0; j<jend; ++j)
1205  result(i, j) = left() * right(i, j);
1206  break;
1207  }
1208  }
1209  }
1210  };
1211  } //namespace
1212 
1213  template<typename DeviceType>
1214  template<typename outputValueType, class ...outputProperties,
1215  typename leftInputValueType, class ...leftInputProperties,
1216  typename rightInputValueType, class ...rightInputProperties>
1217  void
1219  matmatProduct( Kokkos::DynRankView<outputValueType, outputProperties...> output,
1220  const Kokkos::DynRankView<leftInputValueType, leftInputProperties...> leftInput,
1221  const Kokkos::DynRankView<rightInputValueType,rightInputProperties...> rightInput,
1222  const bool hasField,
1223  const bool isTranspose ) {
1224 
1225  typedef Kokkos::DynRankView<outputValueType, outputProperties...> OutputViewType;
1226  typedef const Kokkos::DynRankView<leftInputValueType, leftInputProperties...> leftInputViewType;
1227  typedef const Kokkos::DynRankView<rightInputValueType,rightInputProperties...> rightInputViewType;
1229 
1230  const size_type loopSize = ( hasField ? output.extent(0)*output.extent(1)*output.extent(2) :
1231  output.extent(0)*output.extent(1) );
1232  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
1233  Kokkos::parallel_for( policy, FunctorType(output, leftInput, rightInput, hasField, isTranspose) );
1234  }
1235 
1236 
1237 
1238 
1239  template<typename DeviceType>
1240  template<typename outputFieldValueType, class ...outputFieldProperties,
1241  typename inputDataValueType, class ...inputDataProperties,
1242  typename inputFieldValueType, class ...inputFieldProperties>
1243  void
1245  matmatProductDataField( Kokkos::DynRankView<outputFieldValueType,outputFieldProperties...> outputFields,
1246  const Kokkos::DynRankView<inputDataValueType, inputDataProperties...> inputData,
1247  const Kokkos::DynRankView<inputFieldValueType, inputFieldProperties...> inputFields,
1248  const char transpose ) {
1249 
1250 #ifdef HAVE_INTREPID2_DEBUG
1251  {
1252  /*
1253  * Check array rank and spatial dimension range (if applicable)
1254  * (1) inputData(C,P), (C,P,D) or (C,P,D,D); P=1 is admissible to allow multiply by const. data
1255  * (2) inputFields(C,F,P,D,D) or (F,P,D,D);
1256  * (3) outputFields(C,F,P,D,D)
1257  */
1258  // (1) inputData is (C,P), (C, P, D) or (C, P, D, D) and 1 <= D <= 3 is required
1259  INTREPID2_TEST_FOR_EXCEPTION( inputData.rank() < 2 ||
1260  inputData.rank() > 4, std::invalid_argument,
1261  ">>> ERROR (ArrayTools::matmatProductDataField): inputData must have rank 2,3 or 4" );
1262  if (inputData.rank() > 2) {
1263  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(2) < 1 ||
1264  inputData.extent(2) > 3, std::invalid_argument,
1265  ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension (2) must be 1,2 or 3" );
1266  }
1267  if (inputData.rank() == 4) {
1268  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(3) < 1 ||
1269  inputData.extent(3) > 3, std::invalid_argument,
1270  ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension (3) must be 1,2 or 3" );
1271  }
1272 
1273  // (2) inputFields is (C,F,P,D,D) or (F,P,D,D) and 1 <= (dimension(rank-1), (rank-2)) <= 3 is required.
1274  INTREPID2_TEST_FOR_EXCEPTION( inputFields.rank() < 4 ||
1275  inputFields.rank() > 5, std::invalid_argument,
1276  ">>> ERROR (ArrayTools::matmatProductDataField): inputFields must have rank 4 or 5");
1277  INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(inputFields.rank()-1) < 1 ||
1278  inputFields.extent(inputFields.rank()-1) > 3, std::invalid_argument,
1279  ">>> ERROR (ArrayTools::matmatProductDataField): inputFields dimension (rank-1) must be 1,2 or 3" );
1280  INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(inputFields.rank()-2) < 1 ||
1281  inputFields.extent(inputFields.rank()-2) > 3, std::invalid_argument,
1282  ">>> ERROR (ArrayTools::matmatProductDataField): inputFields dimension (rank-2) must be 1,2 or 3" );
1283 
1284  // (3) outputFields is (C,F,P,D,D)
1285  INTREPID2_TEST_FOR_EXCEPTION( outputFields.rank() != 5, std::invalid_argument,
1286  ">>> ERROR (ArrayTools::matmatProductDataField): outputFields must have rank 5" );
1287  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(3) < 1 ||
1288  outputFields.extent(3) > 3, std::invalid_argument,
1289  ">>> ERROR (ArrayTools::matmatProductDataField): outputFields dimension (3) must be 1,2 or 3" );
1290  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(4) < 1 ||
1291  outputFields.extent(4) > 3, std::invalid_argument,
1292  ">>> ERROR (ArrayTools::matmatProductDataField): outputFields dimension (4) must be 1,2 or 3" );
1293 
1294  /*
1295  * Dimension cross-checks:
1296  * (1) inputData vs. inputFields
1297  * (2) outputFields vs. inputData
1298  * (3) outputFields vs. inputFields
1299  *
1300  * Cross-check (2): outputFields(C,F,P,D,D) vs. inputData(C,P), (C,P,D) or (C,P,D,D):
1301  * dimensions C, and D must match in all cases, dimension P must match only when non-constant
1302  * data is specified (P>1). Do not check P dimensions with constant data, i.e., when P=1 in
1303  * inputData(C,1,...)
1304  */
1305  if (inputData.extent(1) > 1) { // check P dimension if P>1 in inputData
1306  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(2) != inputData.extent(1), std::invalid_argument,
1307  ">>> ERROR (ArrayTools::matmatProductDataField): outputFields dimension (2) does not match to inputData dimension (1)" );
1308  }
1309  if (inputData.rank() == 2) { // inputData(C,P) -> C match
1310  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(0) != inputData.extent(0), std::invalid_argument,
1311  ">>> ERROR (ArrayTools::matmatProductDataField): outputFields dimension (0) does not match to inputData dimension (0)" );
1312  }
1313  if (inputData.rank() == 3) { // inputData(C,P,D) -> C, D match
1314  const size_type f1[] = { 0, 3, 4 }, f2[] = { 0, 2, 2 };
1315  for (size_type i=0; i<3; ++i) {
1316  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputData.extent(f2[i]), std::invalid_argument,
1317  ">>> ERROR (ArrayTools::matmatProductDataField): outputFields dimension does not match to inputData dimension" );
1318  }
1319  }
1320  if (inputData.rank() == 4) { // inputData(C,P,D,D) -> C, D, D match
1321  const size_type f1[] = { 0, 3, 4 }, f2[] = { 0, 2, 3 };
1322  for (size_type i=0; i<3; ++i) {
1323  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputData.extent(f2[i]), std::invalid_argument,
1324  ">>> ERROR (ArrayTools::matmatProductDataField): outputFields dimension does not match to inputData dimension" );
1325  }
1326  }
1327 
1328  /*
1329  * Cross-checks (1,3):
1330  */
1331  if (inputFields.rank() == 5) {
1332  // 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
1333  if (inputData.extent(1) > 1) { // check P dimension if P>1 in inputData
1334  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(1) != inputFields.extent(2), std::invalid_argument,
1335  ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension (1) does not match to inputFields dimension (2)" );
1336  }
1337  if (inputData.rank() == 2) { // inputData(C,P) -> C match
1338  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(0) != inputFields.extent(0), std::invalid_argument,
1339  ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension (0) does not match to inputFields dimension (0)" );
1340  }
1341  if (inputData.rank() == 3) { // inputData(C,P,D) -> C, D match
1342 
1343  const size_type f1[] = { 0, 2, 2 }, f2[] = { 0, 3, 4 };
1344  for (size_type i=0; i<3; ++i) {
1345  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
1346  ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension does not match to inputFields dimension" );
1347  }
1348  }
1349  if (inputData.rank() == 4) { // inputData(C,P,D,D) -> C, D, D match
1350  const size_type f1[] = { 0, 2, 3 }, f2[] = { 0, 3, 4 };
1351  for (size_type i=0; i<3; ++i) {
1352  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
1353  ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension does not match to inputFields dimension" );
1354  }
1355  }
1356 
1357  // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(C,F,P,D,D): all dimensions C, F, P, D must match
1358  for (size_type i=0; i<outputFields.rank(); ++i) {
1359  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(i) != inputFields.extent(i), std::invalid_argument,
1360  ">>> ERROR (ArrayTools::matmatProductDataField): outputFields dimension does not match to inputFields dimension" );
1361  }
1362  } else {
1363  // 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
1364  if (inputData.extent(1) > 1) { // check P if P>1 in inputData
1365  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(1) != inputFields.extent(1), std::invalid_argument,
1366  ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension (1) does not match to inputFields dimension (1)" );
1367  }
1368  if (inputData.rank() == 3) {
1369  const size_type f1[] = { 2, 2 }, f2[] = { 2, 3 };
1370  for (size_type i=0; i<2; ++i) {
1371  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
1372  ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension does not match to inputFields dimension" );
1373  }
1374  }
1375  if (inputData.rank() == 4) {
1376  const size_type f1[] = { 2, 3 }, f2[] = { 2, 3 };
1377  for (size_type i=0; i<2; ++i) {
1378  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
1379  ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension does not match to inputFields dimension" );
1380  }
1381  }
1382 
1383  // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(F,P,D,D): dimensions F, P, D must match
1384  {
1385  const size_type f1[] = { 1, 2, 3, 4 }, f2[] = { 0, 1, 2, 3 };
1386  for (size_type i=0; i<4; ++i) {
1387  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
1388  ">>> ERROR (ArrayTools::matmatProductDataField): outputFields dimension does not match to inputFields dimension" );
1389  }
1390  }
1391  }
1392  }
1393 #endif
1394  // body
1396  inputData,
1397  inputFields,
1398  true,
1399  transpose == 't' || transpose == 'T' );
1400  }
1401 
1402 
1403 
1404  template<typename DeviceType>
1405  template<typename outputDataValueType, class ...outputDataProperties,
1406  typename inputDataLeftValueType, class ...inputDataLeftProperties,
1407  typename inputDataRightValueType, class ...inputDataRightProperties>
1408  void
1409  ArrayTools<DeviceType>::matmatProductDataData( Kokkos::DynRankView<outputDataValueType, outputDataProperties...> outputData,
1410  const Kokkos::DynRankView<inputDataLeftValueType, inputDataLeftProperties...> inputDataLeft,
1411  const Kokkos::DynRankView<inputDataRightValueType,inputDataRightProperties...> inputDataRight,
1412  const char transpose ) {
1413 
1414 #ifdef HAVE_INTREPID2_DEBUG
1415  {
1416  /*
1417  * Check array rank and spatial dimension range (if applicable)
1418  * (1) inputDataLeft(C,P), (C,P,D) or (C,P,D,D); P=1 is admissible to allow multiply by const. left data
1419  * (2) inputDataRight(C,P,D,D) or (P,D,D);
1420  * (3) outputData(C,P,D,D)
1421  */
1422  // (1) inputDataLeft is (C,P), (C,P,D) or (C,P,D,D) and 1 <= D <= 3 is required
1423  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.rank() < 2 ||
1424  inputDataLeft.rank() > 4, std::invalid_argument,
1425  ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft must have rank 2,3 or 4" );
1426  if (inputDataLeft.rank() > 2) {
1427  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(2) < 1 ||
1428  inputDataLeft.extent(2) > 3, std::invalid_argument,
1429  ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft dimension (2) must be 1,2 or 3" );
1430  }
1431  if (inputDataLeft.rank() == 4) {
1432  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(3) < 1 ||
1433  inputDataLeft.extent(3) > 3, std::invalid_argument,
1434  ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft dimension (3) must be 1,2 or 3" );
1435  }
1436 
1437  // (2) inputDataRight is (C,P,D,D) or (P,D,D) and 1 <= (D=dimension(rank-1),(rank-2)) <= 3 is required.
1438  INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.rank() < 3 ||
1439  inputDataRight.rank() > 4, std::invalid_argument,
1440  ">>> ERROR (ArrayTools::matmatProductDataData): inputDataRight must have rank 3 or 4" );
1441  INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.extent(inputDataRight.rank()-1) < 1 ||
1442  inputDataRight.extent(inputDataRight.rank()-1) > 3, std::invalid_argument,
1443  ">>> ERROR (ArrayTools::matmatProductDataData): inputDataRight (rank-1) must be 1,2 or 3" );
1444  INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.extent(inputDataRight.rank()-2) < 1 ||
1445  inputDataRight.extent(inputDataRight.rank()-2) > 3, std::invalid_argument,
1446  ">>> ERROR (ArrayTools::matmatProductDataData): inputDataRight (rank-2) must be 1,2 or 3" );
1447 
1448  // (3) outputData is (C,P,D,D)
1449  INTREPID2_TEST_FOR_EXCEPTION( outputData.rank() != 4, std::invalid_argument,
1450  ">>> ERROR (ArrayTools::matmatProductDataData): outputData must have rank 4" );
1451  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(2) < 1 ||
1452  outputData.extent(2) > 3, std::invalid_argument,
1453  ">>> ERROR (ArrayTools::matmatProductDataData): outputData (2) must be 1,2 or 3" );
1454  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(3) < 1 ||
1455  outputData.extent(3) > 3, std::invalid_argument,
1456  ">>> ERROR (ArrayTools::matmatProductDataData): outputData (3) must be 1,2 or 3" );
1457 
1458  /*
1459  * Dimension cross-checks:
1460  * (1) inputDataLeft vs. inputDataRight
1461  * (2) outputData vs. inputDataLeft
1462  * (3) outputData vs. inputDataRight
1463  *
1464  * Cross-check (2): outputData(C,P,D,D) vs. inputDataLeft(C,P), (C,P,D) or (C,P,D,D):
1465  * dimensions C, and D must match in all cases, dimension P must match only when non-constant
1466  * data is specified (P>1). Do not check P dimensions with constant left data, i.e., when P=1 in
1467  * inputDataLeft(C,1,...)
1468  */
1469  if (inputDataLeft.extent(1) > 1) { // check P dimension if P>1 in inputDataLeft
1470  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(1) != inputDataLeft.extent(1), std::invalid_argument,
1471  ">>> ERROR (ArrayTools::matmatProductDataData): outputData dimension (1) does not match to inputDataLeft dimension (1)" );
1472  }
1473  if (inputDataLeft.rank() == 2) { // inputDataLeft(C,P): check C
1474  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(0) != inputDataLeft.extent(0), std::invalid_argument,
1475  ">>> ERROR (ArrayTools::matmatProductDataData): outputData dimension (0) does not match to inputDataLeft dimension (0)" );
1476  }
1477  if (inputDataLeft.rank() == 3) { // inputDataLeft(C,P,D): check C and D
1478  const size_type f1[] = { 0, 2, 3 }, f2[] = { 0, 2, 2 };
1479  for (size_type i=0; i<3; ++i) {
1480  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataLeft.extent(f2[i]), std::invalid_argument,
1481  ">>> ERROR (ArrayTools::matmatProductDataData): outputData dimension does not match to inputDataLeft dimension" );
1482  }
1483  }
1484  if (inputDataLeft.rank() == 4) { // inputDataLeft(C,P,D,D): check C and D
1485  const size_type f1[] = { 0, 2, 3 }, f2[] = { 0, 2, 3 };
1486  for (size_type i=0; i<3; ++i) {
1487  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataLeft.extent(f2[i]), std::invalid_argument,
1488  ">>> ERROR (ArrayTools::matmatProductDataData): outputData dimension does not match to inputDataLeft dimension" );
1489  }
1490  }
1491 
1492  /*
1493  * Cross-checks (1,3):
1494  */
1495  if (inputDataRight.rank() == 4) {
1496  // 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
1497  if (inputDataLeft.extent(1) > 1) { // check P dimension if P>1 in inputDataLeft
1498  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(1) != inputDataRight.extent(1), std::invalid_argument,
1499  ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft dimension (1) does not match to inputDataRight dimension (1)" );
1500  }
1501  if (inputDataLeft.rank() == 2) { // inputDataLeft(C,P): check C
1502  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(0) != inputDataRight.extent(0), std::invalid_argument,
1503  ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft dimension (0) does not match to inputDataRight dimension (0)" );
1504  }
1505  if (inputDataLeft.rank() == 3) { // inputDataLeft(C,P,D): check C and D
1506  const size_type f1[] = { 0, 2, 2 }, f2[] = { 0, 2, 3 };
1507  for (size_type i=0; i<3; ++i) {
1508  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1509  ">>> ERROR (ArrayTools::matmatProductDataData): intputDataLeft dimension does not match to inputDataRight dimension" );
1510  }
1511  }
1512  if (inputDataLeft.rank() == 4) { // inputDataLeft(C,P,D,D): check C and D
1513  const size_type f1[] = { 0, 2, 3 }, f2[] = { 0, 2, 3 };
1514  for (size_type i=0; i<3; ++i) {
1515  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1516  ">>> ERROR (ArrayTools::matmatProductDataData): intputDataLeft dimension does not match to inputDataRight dimension" );
1517  }
1518  }
1519 
1520  // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(C,P,D,D): all dimensions C, P, D must match
1521  for (size_type i=0; i< outputData.rank(); ++i) {
1522  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(i) != inputDataRight.extent(i), std::invalid_argument,
1523  ">>> ERROR (ArrayTools::matmatProductDataData): outputData dimension does not match to inputDataRight dimension" );
1524  }
1525  } else {
1526  // Cross-check (1): inputDataLeft(C,P), (C,P,D), or (C,P,D,D) vs. inputDataRight(P,D,D): dimensions P, D must match
1527  if (inputDataLeft.extent(1) > 1) { // check P if P>1 in inputData
1528  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(1) != inputDataRight.extent(0), std::invalid_argument,
1529  ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft dimension (1) does not match to inputDataRight dimension (0)" );
1530  }
1531  if (inputDataLeft.rank() == 3) { // inputDataLeft(C,P,D): check D
1532  const size_type f1[] = { 2, 2 }, f2[] = { 1, 2 };
1533  for (size_type i=0; i<2; ++i) {
1534  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1535  ">>> ERROR (ArrayTools::matmatProductDataData): intputDataLeft dimension does not match to inputDataRight dimension" );
1536  }
1537  }
1538  if (inputDataLeft.rank() == 4) { // inputDataLeft(C,P,D,D): check D
1539  const size_type f1[] = { 2, 3 }, f2[] = { 1, 2 };
1540  for (size_type i=0; i<2; ++i) {
1541  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1542  ">>> ERROR (ArrayTools::matmatProductDataData): intputDataLeft dimension does not match to inputDataRight dimension" );
1543  }
1544  }
1545  // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(P,D,D): dimensions P, D must match
1546  {
1547  const size_type f1[] = { 1, 2, 3 }, f2[] = { 0, 1, 2 };
1548  for (size_type i=0; i<3; ++i) {
1549  INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1550  ">>> ERROR (ArrayTools::matmatProductDataData): outputData dimension does not match to inputDataRight dimension" );
1551  }
1552  }
1553  }
1554  }
1555 #endif
1556  // body
1558  inputDataLeft,
1559  inputDataRight,
1560  false,
1561  transpose == 't' || transpose == 'T' );
1562  }
1563 
1564 } // end namespace Intrepid2
1565 #endif
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...
Functor for crossProduct see Intrepid2::ArrayTools for more.
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...
Functor for outerProduct see Intrepid2::ArrayTools for more.
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 ...
Functor for matmatProduct see Intrepid2::ArrayTools for more.
Functor for matvecProduct see Intrepid2::ArrayTools for more.
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...
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...
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...
Utility class that provides methods for higher-order algebraic manipulation of user-defined arrays...
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 ...
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...