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