49 #ifndef __INTREPID2_ARRAYTOOLS_DEF_TENSOR_HPP__
50 #define __INTREPID2_ARRAYTOOLS_DEF_TENSOR_HPP__
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;
65 KOKKOS_INLINE_FUNCTION
67 leftInputViewType leftInput_,
68 rightInputViewType rightInput_,
70 const bool isCrossProd3D_)
71 : _output(output_), _leftInput(leftInput_), _rightInput(rightInput_),
72 _hasField(hasField_), _isCrossProd3D(isCrossProd3D_) {}
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();
80 unrollIndex( cell, field, point,
86 unrollIndex( cell, point,
91 auto result = ( _hasField ? Kokkos::subview(_output, cell, field, point, Kokkos::ALL()) :
92 Kokkos::subview(_output, cell, point, Kokkos::ALL()));
94 auto left = Kokkos::subview(_leftInput, cell, point, Kokkos::ALL());
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()));
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);
108 result(0) = left(0)*right(1) - left(1)*right(0);
114 template<
typename SpT>
115 template<
typename outputValueType,
class ...outputProperties,
116 typename leftInputValueType,
class ...leftInputProperties,
117 typename rightInputValueType,
class ...rightInputProperties>
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 ) {
125 typedef Kokkos::DynRankView<outputValueType, outputProperties...> outputViewType;
126 typedef const Kokkos::DynRankView<leftInputValueType, leftInputProperties...> leftInputViewType;
127 typedef const Kokkos::DynRankView<rightInputValueType,rightInputProperties...> rightInputViewType;
129 typedef typename ExecSpace< typename rightInputViewType::execution_space , SpT >::ExecSpaceType ExecSpaceType;
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) );
140 template<
typename ExecSpaceType>
141 template<
typename outputFieldValueType,
class ...outputFieldProperties,
142 typename inputDataValueType,
class ...inputDataProperties,
143 typename inputFieldValueType,
class ...inputFieldProperties>
147 const Kokkos::DynRankView<inputDataValueType, inputDataProperties...> inputData,
148 const Kokkos::DynRankView<inputFieldValueType, inputFieldProperties...> inputFields ) {
150 #ifdef HAVE_INTREPID2_DEBUG
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");
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" );
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");
182 if (inputFields.rank() == 4) {
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");
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");
200 if (inputData.extent(2) == 2) {
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");
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");
219 if (inputData.extent(2) == 2) {
221 if (inputFields.rank() == 4) {
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");
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");
238 if (inputFields.rank() == 4) {
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");
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");
263 template<
typename ExecSpaceType>
264 template<
typename outputDataValueType,
class ...outputDataProperties,
265 typename inputDataLeftValueType,
class ...inputDataLeftProperties,
266 typename inputDataRightValueType,
class ...inputDataRightProperties>
270 const Kokkos::DynRankView<inputDataLeftValueType, inputDataLeftProperties...> inputDataLeft,
271 const Kokkos::DynRankView<inputDataRightValueType,inputDataRightProperties...> inputDataRight ) {
273 #ifdef HAVE_INTREPID2_DEBUG
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");
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" );
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)");
305 if (inputDataRight.rank() == 3) {
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");
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");
323 if (inputDataLeft.extent(2) == 2) {
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");
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");
340 if (inputDataLeft.extent(2) == 2) {
342 if (inputDataRight.rank() == 3) {
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");
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");
356 if (inputDataRight.rank() == 3) {
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");
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");
381 namespace FunctorArrayTools {
385 template <
typename outputViewType,
typename leftInputViewType,
typename rightInputViewType >
387 outputViewType _output;
388 const leftInputViewType _leftInput;
389 const rightInputViewType _rightInput;
390 const bool _hasField;
392 KOKKOS_INLINE_FUNCTION
394 leftInputViewType leftInput_,
395 rightInputViewType rightInput_,
396 const bool hasField_)
397 : _output(output_), _leftInput(leftInput_), _rightInput(rightInput_),
398 _hasField(hasField_) {}
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();
406 unrollIndex( cell, field, point,
412 unrollIndex( cell, point,
417 auto result = ( _hasField ? Kokkos::subview(_output, cell, field, point, Kokkos::ALL(), Kokkos::ALL()) :
418 Kokkos::subview(_output, cell, point, Kokkos::ALL(), Kokkos::ALL()));
420 auto left = Kokkos::subview(_leftInput, cell, point, Kokkos::ALL());
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()));
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);
437 template<
typename SpT>
438 template<
typename outputValueType,
class ...outputProperties,
439 typename leftInputValueType,
class ...leftInputProperties,
440 typename rightInputValueType,
class ...rightInputProperties>
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 ) {
448 typedef Kokkos::DynRankView<outputValueType, outputProperties...> outputViewType;
449 typedef const Kokkos::DynRankView<leftInputValueType, leftInputProperties...> leftInputViewType;
450 typedef const Kokkos::DynRankView<rightInputValueType,rightInputProperties...> rightInputViewType;
452 typedef typename ExecSpace< typename leftInputViewType::execution_space , SpT >::ExecSpaceType ExecSpaceType;
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) );
461 template<
typename ExecSpaceType>
462 template<
typename outputFieldValueType,
class ...outputFieldProperties,
463 typename inputDataValueType,
class ...inputDataProperties,
464 typename inputFieldValueType,
class ...inputFieldProperties>
468 const Kokkos::DynRankView<inputDataValueType, inputDataProperties...> inputData,
469 const Kokkos::DynRankView<inputFieldValueType, inputFieldProperties...> inputFields ) {
471 #ifdef HAVE_INTREPID2_DEBUG
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");
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" );
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");
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");
521 if (inputFields.rank() == 4) {
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");
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");
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");
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");
568 template<
typename ExecSpaceType>
569 template<
typename outputDataValueType,
class ...outputDataProperties,
570 typename inputDataLeftValuetype,
class ...inputDataLeftProperties,
571 typename inputDataRightValueType,
class ...inputDataRightProperties>
575 const Kokkos::DynRankView<inputDataLeftValuetype, inputDataLeftProperties...> inputDataLeft,
576 const Kokkos::DynRankView<inputDataRightValueType,inputDataRightProperties...> inputDataRight ) {
578 #ifdef HAVE_INTREPID2_DEBUG
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");
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" );
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");
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");
628 if (inputDataRight.rank() == 3) {
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");
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");
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");
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");
672 namespace FunctorArrayTools {
676 template <
typename outputViewType,
677 typename leftInputViewType,
678 typename rightInputViewType>
680 outputViewType _output;
681 const leftInputViewType _leftInput;
682 const rightInputViewType _rightInput;
684 const bool _isTranspose;
686 KOKKOS_INLINE_FUNCTION
688 leftInputViewType leftInput_,
689 rightInputViewType rightInput_,
690 const bool isTranspose_)
691 : _output(output_), _leftInput(leftInput_), _rightInput(rightInput_), _isTranspose(isTranspose_) {}
693 template<
typename resultViewType,
694 typename leftViewType,
695 typename rightViewType>
696 KOKKOS_FORCEINLINE_FUNCTION
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);
705 typedef typename resultViewType::value_type value_type;
707 switch (left.rank()) {
710 for (ordinal_type i=0;i<iend;++i) {
712 for (ordinal_type j=0;j<jend;++j)
713 tmp += left(j, i)*right(j);
717 for (ordinal_type i=0;i<iend;++i) {
719 for (ordinal_type j=0;j<jend;++j)
720 tmp += left(i, j)*right(j);
726 for (ordinal_type i=0;i<iend;++i)
727 result(i) = left(i)*right(i);
731 const value_type val = left();
732 for (ordinal_type i=0;i<iend;++i) {
733 result(i) = val*right(i);
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();
746 auto result = Kokkos::subview(_output, cl, pt, Kokkos::ALL());
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));
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 );
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();
765 auto result = Kokkos::subview(_output, cl, bf, pt, Kokkos::ALL());
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));
772 const auto right = ( rightRank == 3 ? Kokkos::subview(_rightInput, bf, pt, Kokkos::ALL()) :
773 Kokkos::subview(_rightInput, cl, bf, pt, Kokkos::ALL()));
775 apply_matvec_product( result, left, right, _isTranspose );
780 template<
typename SpT>
781 template<
typename outputValueType,
class ...outputProperties,
782 typename leftInputValueType,
class ...leftInputProperties,
783 typename rightInputValueType,
class ...rightInputProperties>
786 matvecProduct( Kokkos::DynRankView<outputValueType, outputProperties...> output,
787 const Kokkos::DynRankView<leftInputValueType, leftInputProperties...> leftInput,
788 const Kokkos::DynRankView<rightInputValueType,rightInputProperties...> rightInput,
790 const bool isTranspose ) {
792 typedef Kokkos::DynRankView<outputValueType, outputProperties...> outputViewType;
793 typedef const Kokkos::DynRankView<leftInputValueType, leftInputProperties...> leftInputViewType;
794 typedef const Kokkos::DynRankView<rightInputValueType,rightInputProperties...> rightInputViewType;
796 typedef typename ExecSpace< typename leftInputViewType::execution_space , SpT >::ExecSpaceType ExecSpaceType;
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) );
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) );
813 template<
typename ExecSpaceType>
814 template<
typename outputFieldValueType,
class ...outputFieldProperties,
815 typename inputDataValueType,
class ...inputDataProperties,
816 typename inputFieldValueType,
class ...inputFieldProperties>
819 const Kokkos::DynRankView<inputDataValueType, inputDataProperties...> inputData,
820 const Kokkos::DynRankView<inputFieldValueType, inputFieldProperties...> inputFields,
821 const char transpose ) {
823 #ifdef HAVE_INTREPID2_DEBUG
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");
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");
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" );
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" );
871 if (inputData.extent(1) > 1) {
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)" );
875 if (inputData.rank() == 2) {
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)" );
879 if (inputData.rank() == 3) {
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" );
886 if (inputData.rank() == 4) {
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" );
897 if (inputFields.rank() == 4) {
899 if (inputData.extent(1) > 1) {
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)" );
903 if (inputData.rank() == 2) {
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)" );
907 if (inputData.rank() == 3) {
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" );
914 if (inputData.rank() == 4) {
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" );
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" );
929 if (inputData.extent(1) > 1) {
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)" );
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)" );
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" );
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" );
961 transpose ==
't' || transpose ==
'T' );
966 template<
typename ExecSpaceType>
967 template<
typename outputDataValueType,
class ...outputDataProperties,
968 typename inputDataLeftValueType,
class ...inputDataLeftProperties,
969 typename inputDataRightValueType,
class ...inputDataRightProperties>
973 const Kokkos::DynRankView<inputDataLeftValueType, inputDataLeftProperties...> inputDataLeft,
974 const Kokkos::DynRankView<inputDataRightValueType,inputDataRightProperties...> inputDataRight,
975 const char transpose ) {
977 #ifdef HAVE_INTREPID2_DEBUG
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" );
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");
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");
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" );
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");
1027 if (inputDataLeft.extent(1) > 1) {
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)" );
1031 if (inputDataLeft.rank() == 2) {
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)" );
1035 if (inputDataLeft.rank() == 3) {
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" );
1042 if (inputDataLeft.rank() == 4) {
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" );
1053 if (inputDataRight.rank() == 3) {
1055 if (inputDataLeft.extent(1) > 1) {
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)" );
1059 if (inputDataLeft.rank() == 2) {
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)" );
1063 if (inputDataLeft.rank() == 3) {
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" );
1070 if (inputDataLeft.rank() == 4) {
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" );
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" );
1085 if (inputDataLeft.extent(1) > 1) {
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)" );
1089 if (inputDataLeft.rank() == 3) {
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)" );
1093 if (inputDataLeft.rank() == 4) {
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" );
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" );
1117 transpose ==
't' || transpose ==
'T' );
1121 namespace FunctorArrayTools {
1125 template <
typename outputViewType,
typename leftInputViewType,
typename rightInputViewType >
1127 outputViewType _output;
1128 leftInputViewType _leftInput;
1129 rightInputViewType _rightInput;
1130 const bool _hasField, _isTranspose;
1131 typedef typename leftInputViewType::value_type value_type;
1133 KOKKOS_INLINE_FUNCTION
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_) {}
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();
1149 unrollIndex( cell, field, point,
1155 unrollIndex( cell, point,
1160 auto result = ( _hasField ? Kokkos::subview(_output, cell, field, point, Kokkos::ALL(), Kokkos::ALL()) :
1161 Kokkos::subview(_output, cell, point, Kokkos::ALL(), Kokkos::ALL()));
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) );
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())));
1175 const ordinal_type iend = result.extent(0);
1176 const ordinal_type jend = result.extent(1);
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);
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);
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);
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);
1216 template<
typename SpT>
1217 template<
typename outputValueType,
class ...outputProperties,
1218 typename leftInputValueType,
class ...leftInputProperties,
1219 typename rightInputValueType,
class ...rightInputProperties>
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 ) {
1228 typedef Kokkos::DynRankView<outputValueType, outputProperties...> outputViewType;
1229 typedef const Kokkos::DynRankView<leftInputValueType, leftInputProperties...> leftInputViewType;
1230 typedef const Kokkos::DynRankView<rightInputValueType,rightInputProperties...> rightInputViewType;
1232 typedef typename ExecSpace< typename leftInputViewType::execution_space , SpT >::ExecSpaceType ExecSpaceType;
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) );
1243 template<
typename ExecSpaceType>
1244 template<
typename outputFieldValueType,
class ...outputFieldProperties,
1245 typename inputDataValueType,
class ...inputDataProperties,
1246 typename inputFieldValueType,
class ...inputFieldProperties>
1250 const Kokkos::DynRankView<inputDataValueType, inputDataProperties...> inputData,
1251 const Kokkos::DynRankView<inputFieldValueType, inputFieldProperties...> inputFields,
1252 const char transpose ) {
1254 #ifdef HAVE_INTREPID2_DEBUG
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" );
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" );
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" );
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" );
1309 if (inputData.extent(1) > 1) {
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)" );
1313 if (inputData.rank() == 2) {
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)" );
1317 if (inputData.rank() == 3) {
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" );
1324 if (inputData.rank() == 4) {
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" );
1335 if (inputFields.rank() == 5) {
1337 if (inputData.extent(1) > 1) {
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)" );
1341 if (inputData.rank() == 2) {
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)" );
1345 if (inputData.rank() == 3) {
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" );
1353 if (inputData.rank() == 4) {
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" );
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" );
1368 if (inputData.extent(1) > 1) {
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)" );
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" );
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" );
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" );
1403 transpose ==
't' || transpose ==
'T' );
1408 template<
typename ExecSpaceType>
1409 template<
typename outputDataValueType,
class ...outputDataProperties,
1410 typename inputDataLeftValueType,
class ...inputDataLeftProperties,
1411 typename inputDataRightValueType,
class ...inputDataRightProperties>
1414 const Kokkos::DynRankView<inputDataLeftValueType, inputDataLeftProperties...> inputDataLeft,
1415 const Kokkos::DynRankView<inputDataRightValueType,inputDataRightProperties...> inputDataRight,
1416 const char transpose ) {
1418 #ifdef HAVE_INTREPID2_DEBUG
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" );
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" );
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" );
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" );
1473 if (inputDataLeft.extent(1) > 1) {
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)" );
1477 if (inputDataLeft.rank() == 2) {
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)" );
1481 if (inputDataLeft.rank() == 3) {
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" );
1488 if (inputDataLeft.rank() == 4) {
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" );
1499 if (inputDataRight.rank() == 4) {
1501 if (inputDataLeft.extent(1) > 1) {
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)" );
1505 if (inputDataLeft.rank() == 2) {
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)" );
1509 if (inputDataLeft.rank() == 3) {
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" );
1516 if (inputDataLeft.rank() == 4) {
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" );
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" );
1531 if (inputDataLeft.extent(1) > 1) {
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)" );
1535 if (inputDataLeft.rank() == 3) {
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" );
1542 if (inputDataLeft.rank() == 4) {
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" );
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" );
1565 transpose ==
't' || transpose ==
'T' );