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 DeviceType>
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;
130 const size_type loopSize = ( hasField ? output.extent(0)*output.extent(1)*output.extent(2) :
131 output.extent(0)*output.extent(1) );
132 const bool isCrossProd3D = (leftInput.extent(2) == 3);
133 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
134 Kokkos::parallel_for( policy, FunctorType(output, leftInput, rightInput, hasField, isCrossProd3D) );
139 template<
typename DeviceType>
140 template<
typename outputFieldValueType,
class ...outputFieldProperties,
141 typename inputDataValueType,
class ...inputDataProperties,
142 typename inputFieldValueType,
class ...inputFieldProperties>
145 crossProductDataField( Kokkos::DynRankView<outputFieldValueType,outputFieldProperties...> outputFields,
146 const Kokkos::DynRankView<inputDataValueType, inputDataProperties...> inputData,
147 const Kokkos::DynRankView<inputFieldValueType, inputFieldProperties...> inputFields ) {
149#ifdef HAVE_INTREPID2_DEBUG
158 INTREPID2_TEST_FOR_EXCEPTION( inputData.rank() != 3, std::invalid_argument,
159 ">>> ERROR (ArrayTools::crossProductDataField): inputData must have rank 3");
160 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(2) < 2 || inputData.extent(2) > 3, std::invalid_argument,
161 ">>> ERROR (ArrayTools::crossProductDataField): inputData dimension(2) must be 2 or 3");
164 INTREPID2_TEST_FOR_EXCEPTION( inputFields.rank() < 3 || inputFields.rank() > 4, std::invalid_argument,
165 ">>> ERROR (ArrayTools::crossProductDataField): inputFields must have rank 3 or 4" );
166 INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(inputFields.rank()-1) < 2 ||
167 inputFields.extent(inputFields.rank()-1) > 3, std::invalid_argument,
168 ">>> ERROR (ArrayTools::crossProductDataField): inputFields dimension (rank-1) must have rank 2 or 3" );
171 INTREPID2_TEST_FOR_EXCEPTION( outputFields.rank() != inputData.extent(2)+1, std::invalid_argument,
172 ">>> ERROR (ArrayTools::crossProductDataField): outputFields rank must match to inputData dimension(2)+1");
181 if (inputFields.rank() == 4) {
183 const size_type f1[] = { 0, 1, 2 }, f2[] = { 0, 2, 3 };
184 for (size_type i=0; i<3; ++i) {
185 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
186 ">>> ERROR (ArrayTools::crossProductDataField): inputData dimension does not match with inputFields");
190 const size_type f1[] = { 1, 2 }, f2[] = { 1, 2 };
191 for (size_type i=0; i<2; ++i) {
192 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
193 ">>> ERROR (ArrayTools::crossProductDataField): inputData dimension does not match with inputFields");
199 if (inputData.extent(2) == 2) {
202 const size_type f1[] = { 0, 2 }, f2[] = { 0, 1 };
203 for (size_type i=0; i<2; ++i) {
204 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputData.extent(f2[i]), std::invalid_argument,
205 ">>> ERROR (ArrayTools::crossProductDataField): outputFields dimension does not match with inputData");
209 const size_type f1[] = { 0, 2, 3 }, f2[] = { 0, 1, 2 };
210 for (size_type i=0; i<3; ++i) {
211 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputData.extent(f2[i]), std::invalid_argument,
212 ">>> ERROR (ArrayTools::crossProductDataField): outputFields dimension does not match with inputData");
218 if (inputData.extent(2) == 2) {
220 if (inputFields.rank() == 4) {
222 const size_type f1[] = { 0, 1, 2 }, f2[] = { 0, 1, 2 };
223 for (size_type i=0; i<3; ++i) {
224 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
225 ">>> ERROR (ArrayTools::crossProductDataField): outputFields dimension does not match with inputFields");
229 const size_type f1[] = { 1, 2 }, f2[] = { 0, 1 };
230 for (size_type i=0; i<2; ++i) {
231 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
232 ">>> ERROR (ArrayTools::crossProductDataField): outputFields dimension does not match with inputFields");
237 if (inputFields.rank() == 4) {
239 for (size_type i=0; i<outputFields.rank(); ++i) {
240 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(i) != inputFields.extent(i), std::invalid_argument,
241 ">>> ERROR (ArrayTools::crossProductDataField): outputFields dimension does not match with inputFields");
245 const size_type f1[] = { 1, 2, 3 }, f2[] = { 0, 1, 2 };
246 for (size_type i=0; i<3; ++i) {
247 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
248 ">>> ERROR (ArrayTools::crossProductDataField): outputFields dimension does not match with inputFields");
262 template<
typename DeviceType>
263 template<
typename outputDataValueType,
class ...outputDataProperties,
264 typename inputDataLeftValueType,
class ...inputDataLeftProperties,
265 typename inputDataRightValueType,
class ...inputDataRightProperties>
268 crossProductDataData( Kokkos::DynRankView<outputDataValueType, outputDataProperties...> outputData,
269 const Kokkos::DynRankView<inputDataLeftValueType, inputDataLeftProperties...> inputDataLeft,
270 const Kokkos::DynRankView<inputDataRightValueType,inputDataRightProperties...> inputDataRight ) {
272#ifdef HAVE_INTREPID2_DEBUG
281 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.rank() != 3, std::invalid_argument,
282 ">>> ERROR (ArrayTools::crossProductDataData): inputDataLeft must have rank 3");
283 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(2) < 2 || inputDataLeft.extent(2) > 3, std::invalid_argument,
284 ">>> ERROR (ArrayTools::crossProductDataData): inputDataLeft dimension(2) must be 2 or 3");
287 INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.rank() < 2 || inputDataRight.rank() > 3, std::invalid_argument,
288 ">>> ERROR (ArrayTools::crossProductDataData): inputDataRight must have rank 2 or 3" );
289 INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.extent(inputDataRight.rank()-1) < 2 ||
290 inputDataRight.extent(inputDataRight.rank()-1) > 3, std::invalid_argument,
291 ">>> ERROR (ArrayTools::crossProductDataData): inputDataRight dimension (rank-1) must have rank 2 or 3" );
294 INTREPID2_TEST_FOR_EXCEPTION( outputData.rank() != inputDataLeft.extent(2), std::invalid_argument,
295 ">>> ERROR (ArrayTools::crossProductDataData): outputData rank must match to inputDataLeft dimension(2)");
304 if (inputDataRight.rank() == 3) {
306 for (size_type i=0; i<inputDataLeft.rank(); ++i) {
307 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(i) != inputDataRight.extent(i), std::invalid_argument,
308 ">>> ERROR (ArrayTools::crossProductDataData): inputDataLeft dimension does not match to inputDataRight");
313 const size_type f1[] = { 1, 2 }, f2[] = { 0, 1 };
314 for (size_type i=0; i<2; ++i) {
315 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
316 ">>> ERROR (ArrayTools::crossProductDataData): inputDataLeft dimension does not match to inputDataRight");
322 if (inputDataLeft.extent(2) == 2) {
324 const size_type f1[] = { 1, 2 }, f2[] = { 0, 1 };
325 for (size_type i=0; i<2; ++i) {
326 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != outputData.extent(f2[i]), std::invalid_argument,
327 ">>> ERROR (ArrayTools::crossProductDataData): inputDataLeft dimension does not match to outputData");
331 for (size_type i=0; i<inputDataLeft.rank(); ++i) {
332 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(i) != outputData.extent(i), std::invalid_argument,
333 ">>> ERROR (ArrayTools::crossProductDataData): inputDataLeft dimension does not match to outputData");
339 if (inputDataLeft.extent(2) == 2) {
341 if (inputDataRight.rank() == 3) {
343 const size_type f1[] = { 0, 1 }, f2[] = { 0, 1 };
344 for (size_type i=0; i<2; ++i) {
345 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
346 ">>> ERROR (ArrayTools::crossProductDataData): outputData dimension does not match to inputDataRight");
350 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(1) != inputDataRight.extent(0), std::invalid_argument,
351 ">>> ERROR (ArrayTools::crossProductDataData): outputData dimension does not match to inputDataRight");
355 if (inputDataRight.rank() == 3) {
357 for (size_type i=0; i<outputData.rank(); ++i) {
358 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(i) != inputDataRight.extent(i), std::invalid_argument,
359 ">>> ERROR (ArrayTools::crossProductDataData): outputData dimension does not match to inputDataRight");
363 const size_type f1[] = { 1, 2 }, f2[] = { 0, 1 };
364 for (size_type i=0; i<2; ++i) {
365 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
366 ">>> ERROR (ArrayTools::crossProductDataData): outputData dimension does not match to inputDataRight");
380 namespace FunctorArrayTools {
384 template <
typename OutputViewType,
typename leftInputViewType,
typename rightInputViewType >
386 OutputViewType _output;
387 const leftInputViewType _leftInput;
388 const rightInputViewType _rightInput;
389 const bool _hasField;
391 KOKKOS_INLINE_FUNCTION
393 leftInputViewType leftInput_,
394 rightInputViewType rightInput_,
395 const bool hasField_)
396 : _output(output_), _leftInput(leftInput_), _rightInput(rightInput_),
397 _hasField(hasField_) {}
399 KOKKOS_INLINE_FUNCTION
400 void operator()(
const size_type iter)
const {
401 size_type cell, field(0), point;
402 size_type rightRank = _rightInput.rank();
405 unrollIndex( cell, field, point,
411 unrollIndex( cell, point,
416 auto result = ( _hasField ? Kokkos::subview(_output, cell, field, point, Kokkos::ALL(), Kokkos::ALL()) :
417 Kokkos::subview(_output, cell, point, Kokkos::ALL(), Kokkos::ALL()));
419 auto left = Kokkos::subview(_leftInput, cell, point, Kokkos::ALL());
421 auto right = (rightRank == 2 + size_type(_hasField)) ?
422 ( _hasField ? Kokkos::subview(_rightInput, field, point, Kokkos::ALL()) :
423 Kokkos::subview(_rightInput, point, Kokkos::ALL())) :
424 ( _hasField ? Kokkos::subview(_rightInput, cell, field, point, Kokkos::ALL()) :
425 Kokkos::subview(_rightInput, cell, point, Kokkos::ALL()));
427 const ordinal_type iend = result.extent(0);
428 const ordinal_type jend = result.extent(1);
429 for (ordinal_type i=0; i<iend; ++i)
430 for (ordinal_type j=0; j<jend; ++j)
431 result(i, j) = left(i)*right(j);
436 template<
typename DeviceType>
437 template<
typename outputValueType,
class ...outputProperties,
438 typename leftInputValueType,
class ...leftInputProperties,
439 typename rightInputValueType,
class ...rightInputProperties>
442 outerProduct( Kokkos::DynRankView<outputValueType, outputProperties...> output,
443 const Kokkos::DynRankView<leftInputValueType, leftInputProperties...> leftInput,
444 const Kokkos::DynRankView<rightInputValueType,rightInputProperties...> rightInput,
445 const bool hasField ) {
447 typedef Kokkos::DynRankView<outputValueType, outputProperties...> OutputViewType;
448 typedef const Kokkos::DynRankView<leftInputValueType, leftInputProperties...> leftInputViewType;
449 typedef const Kokkos::DynRankView<rightInputValueType,rightInputProperties...> rightInputViewType;
452 const size_type loopSize = ( hasField ? output.extent(0)*output.extent(1)*output.extent(2) :
453 output.extent(0)*output.extent(1) );
454 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
455 Kokkos::parallel_for( policy, FunctorType(output, leftInput, rightInput, hasField) );
459 template<
typename DeviceType>
460 template<
typename outputFieldValueType,
class ...outputFieldProperties,
461 typename inputDataValueType,
class ...inputDataProperties,
462 typename inputFieldValueType,
class ...inputFieldProperties>
465 outerProductDataField( Kokkos::DynRankView<outputFieldValueType,outputFieldProperties...> outputFields,
466 const Kokkos::DynRankView<inputDataValueType, inputDataProperties...> inputData,
467 const Kokkos::DynRankView<inputFieldValueType, inputFieldProperties...> inputFields ) {
469#ifdef HAVE_INTREPID2_DEBUG
478 INTREPID2_TEST_FOR_EXCEPTION( inputData.rank() != 3, std::invalid_argument,
479 ">>> ERROR (ArrayTools::outerProductDataField): inputData must have rank 3");
480 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(2) < 2 || inputData.extent(2) > 3, std::invalid_argument,
481 ">>> ERROR (ArrayTools::outerProductDataField): inputData dimension(2) must be 2 or 3");
484 INTREPID2_TEST_FOR_EXCEPTION( inputFields.rank() < 3 || inputFields.rank() > 4, std::invalid_argument,
485 ">>> ERROR (ArrayTools::outerProductDataField): inputFields must have rank 3 or 4" );
486 INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(inputFields.rank()-1) < 2 ||
487 inputFields.extent(inputFields.rank()-1) < 3, std::invalid_argument,
488 ">>> ERROR (ArrayTools::outerProductDataField): inputFields dimension (rank-1) must have rank 2 or 3" );
491 INTREPID2_TEST_FOR_EXCEPTION( outputFields.rank() != 5, std::invalid_argument,
492 ">>> ERROR (ArrayTools::outerProductDataField): outputFields must have rank 5");
493 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(3) < 2 ||
494 outputFields.extent(3) > 3, std::invalid_argument,
495 ">>> ERROR (ArrayTools::outerProductDataField): outputFields dimension(3) must be 2 or 3");
496 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(4) < 2 ||
497 outputFields.extent(4) > 3, std::invalid_argument,
498 ">>> ERROR (ArrayTools::outerProductDataField): outputFields dimension(4) must be 2 or 3");
509 const size_type f1[] = { 0, 2, 3, 4 }, f2[] = { 0, 1, 2, 2 };
510 for (size_type i=0; i<4; ++i) {
511 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputData.extent(f2[i]), std::invalid_argument,
512 ">>> ERROR (ArrayTools::outerProductDataField): outputFields dimension does not match with inputData");
519 if (inputFields.rank() == 4) {
522 const size_type f1[] = { 0, 1, 2 }, f2[] = { 0, 2, 3 };
523 for (size_type i=0; i<3; ++i) {
524 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
525 ">>> ERROR (ArrayTools::outerProductDataField): inputData dimension does not match with inputFields");
531 const size_type f1[] = { 0, 1, 2, 3, 4 }, f2[] = { 0, 1, 2, 3, 3 };
532 for (size_type i=0; i<5; ++i) {
533 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
534 ">>> ERROR (ArrayTools::outerProductDataField): outputFields dimension does not match with inputFields");
540 const size_type f1[] = { 1, 2 }, f2[] = { 1, 2 };
541 for (size_type i=0; i<2; ++i) {
542 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
543 ">>> ERROR (ArrayTools::outerProductDataField): inputData dimension does not match with inputFields");
549 const size_type f1[] = { 1, 2, 3, 4 }, f2[] = { 0, 1, 2, 2 };
550 for (size_type i=0; i<4; ++i) {
551 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
552 ">>> ERROR (ArrayTools::outerProductDataField): outputFields dimension does not match with inputFields");
566 template<
typename DeviceType>
567 template<
typename outputDataValueType,
class ...outputDataProperties,
568 typename inputDataLeftValuetype,
class ...inputDataLeftProperties,
569 typename inputDataRightValueType,
class ...inputDataRightProperties>
572 outerProductDataData( Kokkos::DynRankView<outputDataValueType, outputDataProperties...> outputData,
573 const Kokkos::DynRankView<inputDataLeftValuetype, inputDataLeftProperties...> inputDataLeft,
574 const Kokkos::DynRankView<inputDataRightValueType,inputDataRightProperties...> inputDataRight ) {
576#ifdef HAVE_INTREPID2_DEBUG
585 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.rank() != 3, std::invalid_argument,
586 ">>> ERROR (ArrayTools::outerProductDataField): inputDataLeft must have rank 3");
587 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(2) < 2 || inputDataLeft.extent(2) > 3, std::invalid_argument,
588 ">>> ERROR (ArrayTools::outerProductDataField): inputDataLeft dimension(2) must be 2 or 3");
591 INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.rank() < 2 || inputDataRight.rank() > 3, std::invalid_argument,
592 ">>> ERROR (ArrayTools::outerProductDataField): inputDataRight must have rank 2 or 3" );
593 INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.extent(inputDataRight.rank()-1) < 2 ||
594 inputDataRight.extent(inputDataRight.rank()-1) < 3, std::invalid_argument,
595 ">>> ERROR (ArrayTools::outerProductDataField): inputDataRight dimension (rank-1) must have rank 2 or 3" );
598 INTREPID2_TEST_FOR_EXCEPTION( outputData.rank() != 4, std::invalid_argument,
599 ">>> ERROR (ArrayTools::outerProductDataField): outputData must have rank 5");
600 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(2) < 2 ||
601 outputData.extent(2) > 3, std::invalid_argument,
602 ">>> ERROR (ArrayTools::outerProductDataField): outputData dimension(3) must be 2 or 3");
603 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(3) < 2 ||
604 outputData.extent(3) > 3, std::invalid_argument,
605 ">>> ERROR (ArrayTools::outerProductDataField): outputData dimension(4) must be 2 or 3");
616 const size_type f1[] = { 0, 1, 2, 3 }, f2[] = { 0, 1, 2, 2 };
617 for (size_type i=0; i<4; ++i) {
618 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataLeft.extent(f2[i]), std::invalid_argument,
619 ">>> ERROR (ArrayTools::outerProductDataField): outputData dimension does not match with inputDataLeft");
626 if (inputDataRight.rank() == 3) {
628 for (size_type i=0; i<inputDataLeft.rank(); ++i) {
629 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(i) != inputDataRight.extent(i), std::invalid_argument,
630 ">>> ERROR (ArrayTools::outerProductDataField): inputDataLeft dimension does not match with inputDataRight");
635 const size_type f1[] = { 0, 1, 2, 3 }, f2[] = { 0, 1, 2, 2 };
636 for (size_type i=0; i<4; ++i) {
637 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
638 ">>> ERROR (ArrayTools::outerProductDataField): outputData dimension does not match with inputDataRight");
644 const size_type f1[] = { 1, 2 }, f2[] = { 0, 1 };
645 for (size_type i=0; i<2; ++i) {
646 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
647 ">>> ERROR (ArrayTools::outerProductDataField): inputDataLeft dimension does not match with inputDataRight");
653 const size_type f1[] = { 1, 2, 3 }, f2[] = { 0, 1, 1 };
654 for (size_type i=0; i<3; ++i) {
655 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
656 ">>> ERROR (ArrayTools::outerProductDataField): outputData dimension does not match with inputDataRight");
670 namespace FunctorArrayTools {
674 template <
typename OutputViewType,
675 typename leftInputViewType,
676 typename rightInputViewType>
678 OutputViewType _output;
679 const leftInputViewType _leftInput;
680 const rightInputViewType _rightInput;
682 const bool _isTranspose;
684 KOKKOS_INLINE_FUNCTION
686 leftInputViewType leftInput_,
687 rightInputViewType rightInput_,
688 const bool isTranspose_)
689 : _output(output_), _leftInput(leftInput_), _rightInput(rightInput_), _isTranspose(isTranspose_) {}
691 template<
typename resultViewType,
692 typename leftViewType,
693 typename rightViewType>
694 KOKKOS_FORCEINLINE_FUNCTION
696 apply_matvec_product( resultViewType &result,
697 const leftViewType &left,
698 const rightViewType &right,
699 const bool isTranspose) {
700 const ordinal_type iend = result.extent(0);
701 const ordinal_type jend = right.extent(0);
703 typedef typename resultViewType::value_type value_type;
705 switch (left.rank()) {
708 for (ordinal_type i=0;i<iend;++i) {
710 for (ordinal_type j=0;j<jend;++j)
711 tmp += left(j, i)*right(j);
715 for (ordinal_type i=0;i<iend;++i) {
717 for (ordinal_type j=0;j<jend;++j)
718 tmp += left(i, j)*right(j);
724 for (ordinal_type i=0;i<iend;++i)
725 result(i) = left(i)*right(i);
729 const value_type val = left();
730 for (ordinal_type i=0;i<iend;++i) {
731 result(i) = val*right(i);
738 KOKKOS_INLINE_FUNCTION
739 void operator()(
const ordinal_type cl,
740 const ordinal_type pt)
const {
741 const auto rightRank = _rightInput.rank();
742 const auto leftRank = _leftInput.rank();
744 auto result = Kokkos::subview(_output, cl, pt, Kokkos::ALL());
746 const auto lpt = (_leftInput.extent(1) == 1 ? size_type(0) : pt);
747 const auto left = ( leftRank == 4 ? Kokkos::subview(_leftInput, cl, lpt, Kokkos::ALL(), Kokkos::ALL()) :
748 leftRank == 3 ? Kokkos::subview(_leftInput, cl, lpt, Kokkos::ALL()) :
749 Kokkos::subview(_leftInput, cl, lpt));
751 const auto right = ( rightRank == 2 ? Kokkos::subview(_rightInput, pt, Kokkos::ALL()) :
752 Kokkos::subview(_rightInput, cl, pt, Kokkos::ALL()) );
753 apply_matvec_product( result, left, right, _isTranspose );
756 KOKKOS_INLINE_FUNCTION
757 void operator()(
const ordinal_type cl,
758 const ordinal_type bf,
759 const ordinal_type pt)
const {
760 const auto rightRank = _rightInput.rank();
761 const auto leftRank = _leftInput.rank();
763 auto result = Kokkos::subview(_output, cl, bf, pt, Kokkos::ALL());
765 const auto lpt = (_leftInput.extent(1) == 1 ? size_type(0) : pt);
766 const auto left = ( leftRank == 4 ? Kokkos::subview(_leftInput, cl, lpt, Kokkos::ALL(), Kokkos::ALL()) :
767 leftRank == 3 ? Kokkos::subview(_leftInput, cl, lpt, Kokkos::ALL()) :
768 Kokkos::subview(_leftInput, cl, lpt));
770 const auto right = ( rightRank == 3 ? Kokkos::subview(_rightInput, bf, pt, Kokkos::ALL()) :
771 Kokkos::subview(_rightInput, cl, bf, pt, Kokkos::ALL()));
773 apply_matvec_product( result, left, right, _isTranspose );
778 template<
typename DeviceType>
779 template<
typename outputValueType,
class ...outputProperties,
780 typename leftInputValueType,
class ...leftInputProperties,
781 typename rightInputValueType,
class ...rightInputProperties>
784 matvecProduct( Kokkos::DynRankView<outputValueType, outputProperties...> output,
785 const Kokkos::DynRankView<leftInputValueType, leftInputProperties...> leftInput,
786 const Kokkos::DynRankView<rightInputValueType,rightInputProperties...> rightInput,
788 const bool isTranspose ) {
790 typedef Kokkos::DynRankView<outputValueType, outputProperties...> OutputViewType;
791 typedef const Kokkos::DynRankView<leftInputValueType, leftInputProperties...> leftInputViewType;
792 typedef const Kokkos::DynRankView<rightInputValueType,rightInputProperties...> rightInputViewType;
796 using range_policy_type = Kokkos::MDRangePolicy
797 < ExecSpaceType, Kokkos::Rank<3>, Kokkos::IndexType<ordinal_type> >;
798 range_policy_type policy( { 0, 0, 0 },
799 { output.extent(0), output.extent(1), output.extent(2) } );
800 Kokkos::parallel_for( policy, FunctorType(output, leftInput, rightInput, isTranspose) );
802 using range_policy_type = Kokkos::MDRangePolicy
803 < ExecSpaceType, Kokkos::Rank<2>, Kokkos::IndexType<ordinal_type> >;
804 range_policy_type policy( { 0, 0 },
805 { output.extent(0), output.extent(1) } );
806 Kokkos::parallel_for( policy, FunctorType(output, leftInput, rightInput, isTranspose) );
810 template<
typename DeviceType>
811 template<
typename outputFieldValueType,
class ...outputFieldProperties,
812 typename inputDataValueType,
class ...inputDataProperties,
813 typename inputFieldValueType,
class ...inputFieldProperties>
816 const Kokkos::DynRankView<inputDataValueType, inputDataProperties...> inputData,
817 const Kokkos::DynRankView<inputFieldValueType, inputFieldProperties...> inputFields,
818 const char transpose ) {
820#ifdef HAVE_INTREPID2_DEBUG
829 INTREPID2_TEST_FOR_EXCEPTION( inputData.rank() < 2 || inputData.rank() > 4, std::invalid_argument,
830 ">>> ERROR (ArrayTools::matvecProductDataField): inputData must have rank 2,3 or 4" );
831 if (inputData.rank() > 2) {
832 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(2) < 1 ||
833 inputData.extent(2) > 3, std::invalid_argument,
834 ">>> ERROR (ArrayTools::matvecProductDataField): inputData dimension(2) must be 1,2 or 3");
836 if (inputData.rank() == 4) {
837 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(3) < 1 ||
838 inputData.extent(3) > 3, std::invalid_argument,
839 ">>> ERROR (ArrayTools::matvecProductDataField): inputData dimension(3) must be 1,2 or 3");
843 INTREPID2_TEST_FOR_EXCEPTION( inputFields.rank() < 3 ||
844 inputFields.rank() > 4, std::invalid_argument,
845 ">>> ERROR (ArrayTools::matvecProductDataField): inputFields must have rank 3 or 4" );
846 INTREPID2_TEST_FOR_EXCEPTION( (inputFields.rank()-1) < 1 ||
847 (inputFields.rank()-1) > 3, std::invalid_argument,
848 ">>> ERROR (ArrayTools::matvecProductDataField): inputFields dimension(rank-1) must be 1,2, or 3" );
851 INTREPID2_TEST_FOR_EXCEPTION( outputFields.rank() != 4, std::invalid_argument,
852 ">>> ERROR (ArrayTools::matvecProductDataField): outputFields must have rank 4" );
853 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(3) < 1 ||
854 outputFields.extent(3) > 3, std::invalid_argument,
855 ">>> ERROR (ArrayTools::matvecProductDataField): outputFields dimension(3) must be 1,2 or 3" );
868 if (inputData.extent(1) > 1) {
869 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(2) != inputData.extent(1), std::invalid_argument,
870 ">>> ERROR (ArrayTools::matvecProductDataField): outputFields dimension(2) must match to inputData dimension(1)" );
872 if (inputData.rank() == 2) {
873 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(0) != inputData.extent(0), std::invalid_argument,
874 ">>> ERROR (ArrayTools::matvecProductDataField): outputFields dimension(0) must match to inputData dimension(0)" );
876 if (inputData.rank() == 3) {
877 const size_type f1[] = { 0, 3 }, f2[] = { 0, 2 };
878 for (size_type i=0; i<2; ++i) {
879 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputData.extent(f2[i]), std::invalid_argument,
880 ">>> ERROR (ArrayTools::matvecProductDataField): outputFields dimension does not match to inputData dimension" );
883 if (inputData.rank() == 4) {
884 const size_type f1[] = { 0, 3, 3 }, f2[] = { 0, 2, 3 };
885 for (size_type i=0; i<3; ++i) {
886 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputData.extent(f2[i]), std::invalid_argument,
887 ">>> ERROR (ArrayTools::matvecProductDataField): outputFields dimension does not match to inputData dimension" );
894 if (inputFields.rank() == 4) {
896 if (inputData.extent(1) > 1) {
897 INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(2) != inputData.extent(1), std::invalid_argument,
898 ">>> ERROR (ArrayTools::matvecProductDataField): inputFields dimension (2) does not match to inputData dimension (1)" );
900 if (inputData.rank() == 2) {
901 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(0) != inputFields.extent(0), std::invalid_argument,
902 ">>> ERROR (ArrayTools::matvecProductDataField): inputData dimension (0) does not match to inputFields dimension (0)" );
904 if (inputData.rank() == 3) {
905 const size_type f1[] = { 0, 2 }, f2[] = { 0, 3 };
906 for (size_type i=0; i<2; ++i) {
907 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
908 ">>> ERROR (ArrayTools::matvecProductDataField): inputData dimension does not match to inputFields dimension" );
911 if (inputData.rank() == 4) {
912 const size_type f1[] = { 0, 2, 3 }, f2[] = { 0, 3, 3 };
913 for (size_type i=0; i<3; ++i) {
914 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
915 ">>> ERROR (ArrayTools::matvecProductDataField): inputData dimension does not match to inputFields dimension" );
920 for (size_type i=0; i<outputFields.rank(); ++i) {
921 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(i) != inputFields.extent(i), std::invalid_argument,
922 ">>> ERROR (ArrayTools::matvecProductDataField): outputFields dimension does not match to inputFields dimension" );
926 if (inputData.extent(1) > 1) {
927 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(1) != inputFields.extent(1), std::invalid_argument,
928 ">>> ERROR (ArrayTools::matvecProductDataField): inputData dimension(1) does not match to inputFields dimension (1)" );
930 if (inputData.rank() == 3) {
931 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(2) != inputFields.extent(2), std::invalid_argument,
932 ">>> ERROR (ArrayTools::matvecProductDataField): inputData dimension(2) does not match to inputFields dimension (2)" );
934 if (inputData.rank() == 4) {
935 const size_type f1[] = { 2, 3 }, f2[] = { 2, 2 };
936 for (size_type i=0; i<2; ++i) {
937 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
938 ">>> ERROR (ArrayTools::matvecProductDataField): inputData dimension does not match to inputFields dimension" );
944 const size_type f1[] = { 1, 2, 3 }, f2[] = { 0, 1, 2 };
945 for (size_type i=0; i<3; ++i) {
946 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
947 ">>> ERROR (ArrayTools::matvecProductDataField): outputFields dimension does not match to inputFields dimension" );
958 transpose ==
't' || transpose ==
'T' );
963 template<
typename DeviceType>
964 template<
typename outputDataValueType,
class ...outputDataProperties,
965 typename inputDataLeftValueType,
class ...inputDataLeftProperties,
966 typename inputDataRightValueType,
class ...inputDataRightProperties>
970 const Kokkos::DynRankView<inputDataLeftValueType, inputDataLeftProperties...> inputDataLeft,
971 const Kokkos::DynRankView<inputDataRightValueType,inputDataRightProperties...> inputDataRight,
972 const char transpose ) {
974#ifdef HAVE_INTREPID2_DEBUG
983 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.rank() < 2 ||
984 inputDataLeft.rank() > 4, std::invalid_argument,
985 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft must have rank 2,3 or 4" );
987 if (inputDataLeft.rank() > 2) {
988 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(2) < 1 ||
989 inputDataLeft.extent(2) > 3, std::invalid_argument,
990 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension(2) must be 1, 2 or 3");
992 if (inputDataLeft.rank() == 4) {
993 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(3) < 1 ||
994 inputDataLeft.extent(3) > 3, std::invalid_argument,
995 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension(3) must be 1, 2 or 3");
999 INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.rank() < 2 ||
1000 inputDataRight.rank() > 3, std::invalid_argument,
1001 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataRight must have rank 2 or 3" );
1002 INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.extent(inputDataRight.rank()-1) < 1 ||
1003 inputDataRight.extent(inputDataRight.rank()-1) > 3, std::invalid_argument,
1004 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataRight dimension (rank-1) must be 1,2 or 3" );
1007 INTREPID2_TEST_FOR_EXCEPTION( outputData.rank() != 3, std::invalid_argument,
1008 ">>> ERROR (ArrayTools::matvecProductDataData): outputData must have rank 3" );
1009 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(2) < 1 ||
1010 outputData.extent(2) > 3, std::invalid_argument,
1011 ">>> ERROR (ArrayTools::matvecProductDataData): outputData dimension(2) must be 1, 2 or 3");
1024 if (inputDataLeft.extent(1) > 1) {
1025 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(1) != inputDataLeft.extent(1), std::invalid_argument,
1026 ">>> ERROR (ArrayTools::matvecProductDataData): outputData dimension(1) muat match inputDataLeft dimension(1)" );
1028 if (inputDataLeft.rank() == 2) {
1029 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(0) != inputDataLeft.extent(0), std::invalid_argument,
1030 ">>> ERROR (ArrayTools::matvecProductDataData): outputData dimension(0) muat match inputDataLeft dimension(0)" );
1032 if (inputDataLeft.rank() == 3) {
1033 const size_type f1[] = { 0, 2 }, f2[] = { 0, 2 };
1034 for (size_type i=0; i<2; ++i) {
1035 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataLeft.extent(f2[i]), std::invalid_argument,
1036 ">>> ERROR (ArrayTools::matvecProductDataData): outputData dimension muat match inputDataLeft dimension" );
1039 if (inputDataLeft.rank() == 4) {
1040 size_type f1[] = { 0, 2}, f2[] = { 0, 2};
1041 if (transpose) f2[1] = 3;
1042 for (size_type i=0; i<2; ++i) {
1043 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataLeft.extent(f2[i]), std::invalid_argument,
1044 ">>> ERROR (ArrayTools::matvecProductDataData): outputData dimension muat match inputDataLeft dimension" );
1051 if (inputDataRight.rank() == 3) {
1053 if (inputDataLeft.extent(1) > 1) {
1054 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(1) != inputDataRight.extent(1), std::invalid_argument,
1055 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension (1) muat match to inputDataRight dimension (1)" );
1057 if (inputDataLeft.rank() == 2) {
1058 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(0) != inputDataRight.extent(0), std::invalid_argument,
1059 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension (0) muat match to inputDataRight dimension (0)" );
1061 if (inputDataLeft.rank() == 3) {
1062 const size_type f1[] = { 0, 2 }, f2[] = { 0, 2 };
1063 for (size_type i=0; i<2; ++i) {
1064 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1065 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension muat match to inputDataRight dimension" );
1068 if (inputDataLeft.rank() == 4) {
1069 size_type f1[] = { 0, 3}, f2[] = { 0, 2};
1070 if (transpose) f1[1] = 2;
1071 for (size_type i=0; i<2; ++i) {
1072 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1073 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension muat match to inputDataRight dimension" );
1078 for (size_type i=0; i<outputData.rank()-1; ++i) {
1079 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(i) != inputDataRight.extent(i), std::invalid_argument,
1080 ">>> ERROR (ArrayTools::matvecProductDataData): outputData dimension muat match to inputDataRight dimension" );
1084 if (inputDataLeft.extent(1) > 1) {
1085 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(1) != inputDataRight.extent(0), std::invalid_argument,
1086 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension (1) does mot match to inputDataright dimension (0)" );
1088 if (inputDataLeft.rank() == 3) {
1089 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(2) != inputDataRight.extent(1), std::invalid_argument,
1090 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension (2) does mot match to inputDataright dimension (1)" );
1092 if (inputDataLeft.rank() == 4) {
1093 const size_type f1[] = { 2, 3 }, f2[] = { 1, 1 };
1094 for (size_type i=0; i<2; ++i) {
1095 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1096 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension muat match to inputDataRight dimension" );
1102 const size_type f1[] = { 1, 2 }, f2[] = { 0, 1 };
1103 for (size_type i=0; i<2; ++i) {
1104 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1105 ">>> ERROR (ArrayTools::matvecProductDataData): outputData dimension muat match to inputDataRight dimension" );
1116 transpose ==
't' || transpose ==
'T' );
1120 namespace FunctorArrayTools {
1124 template <
typename OutputViewType,
typename leftInputViewType,
typename rightInputViewType >
1126 OutputViewType _output;
1127 leftInputViewType _leftInput;
1128 rightInputViewType _rightInput;
1129 const bool _hasField, _isTranspose;
1130 typedef typename leftInputViewType::value_type value_type;
1132 KOKKOS_INLINE_FUNCTION
1134 leftInputViewType leftInput_,
1135 rightInputViewType rightInput_,
1136 const bool hasField_,
1137 const bool isTranspose_)
1138 : _output(output_), _leftInput(leftInput_), _rightInput(rightInput_),
1139 _hasField(hasField_), _isTranspose(isTranspose_) {}
1141 KOKKOS_INLINE_FUNCTION
1142 void operator()(
const size_type iter)
const {
1143 size_type cell(0), field(0), point(0);
1144 size_type leftRank = _leftInput.rank();
1145 size_type rightRank = _rightInput.rank();
1148 unrollIndex( cell, field, point,
1154 unrollIndex( cell, point,
1159 auto result = ( _hasField ? Kokkos::subview(_output, cell, field, point, Kokkos::ALL(), Kokkos::ALL()) :
1160 Kokkos::subview(_output, cell, point, Kokkos::ALL(), Kokkos::ALL()));
1162 const auto lpoint = (_leftInput.extent(1) == 1 ? 0 : point);
1163 auto left = ( leftRank == 4 ? Kokkos::subview(_leftInput, cell, lpoint, Kokkos::ALL(), Kokkos::ALL()) :
1164 leftRank == 3 ? Kokkos::subview(_leftInput, cell, lpoint, Kokkos::ALL()) :
1165 Kokkos::subview(_leftInput, cell, lpoint) );
1168 const bool hasCell = (_hasField ? rightRank == 5 : rightRank == 4);
1169 auto right = ( _hasField ? ( hasCell ? Kokkos::subview(_rightInput, cell, field, point, Kokkos::ALL(), Kokkos::ALL()) :
1170 Kokkos::subview(_rightInput, field, point, Kokkos::ALL(), Kokkos::ALL())):
1171 ( hasCell ? Kokkos::subview(_rightInput, cell, point, Kokkos::ALL(), Kokkos::ALL()) :
1172 Kokkos::subview(_rightInput, point, Kokkos::ALL(), Kokkos::ALL())));
1174 const ordinal_type iend = result.extent(0);
1175 const ordinal_type jend = result.extent(1);
1180 const size_type kend = right.extent(0);
1181 for (ordinal_type i=0; i<iend; ++i)
1182 for (ordinal_type j=0; j<jend; ++j) {
1183 result(i, j) = value_type(0);
1184 for (size_type k=0; k<kend; ++k)
1185 result(i, j) += left(k, i) * right(k, j);
1188 const size_type kend = right.extent(0);
1189 for (ordinal_type i=0; i<iend; ++i)
1190 for (ordinal_type j=0; j<jend; ++j) {
1191 result(i, j) = value_type(0);
1192 for (size_type k=0; k<kend; ++k)
1193 result(i, j) += left(i, k) * right(k, j);
1199 for (ordinal_type i=0; i<iend; ++i)
1200 for (ordinal_type j=0; j<jend; ++j)
1201 result(i, j) = left(i) * right(i, j);
1205 for (ordinal_type i=0; i<iend; ++i)
1206 for (ordinal_type j=0; j<jend; ++j)
1207 result(i, j) = left() * right(i, j);
1215 template<
typename DeviceType>
1216 template<
typename outputValueType,
class ...outputProperties,
1217 typename leftInputValueType,
class ...leftInputProperties,
1218 typename rightInputValueType,
class ...rightInputProperties>
1221 matmatProduct( Kokkos::DynRankView<outputValueType, outputProperties...> output,
1222 const Kokkos::DynRankView<leftInputValueType, leftInputProperties...> leftInput,
1223 const Kokkos::DynRankView<rightInputValueType,rightInputProperties...> rightInput,
1224 const bool hasField,
1225 const bool isTranspose ) {
1227 typedef Kokkos::DynRankView<outputValueType, outputProperties...> OutputViewType;
1228 typedef const Kokkos::DynRankView<leftInputValueType, leftInputProperties...> leftInputViewType;
1229 typedef const Kokkos::DynRankView<rightInputValueType,rightInputProperties...> rightInputViewType;
1232 const size_type loopSize = ( hasField ? output.extent(0)*output.extent(1)*output.extent(2) :
1233 output.extent(0)*output.extent(1) );
1234 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
1235 Kokkos::parallel_for( policy, FunctorType(output, leftInput, rightInput, hasField, isTranspose) );
1241 template<
typename DeviceType>
1242 template<
typename outputFieldValueType,
class ...outputFieldProperties,
1243 typename inputDataValueType,
class ...inputDataProperties,
1244 typename inputFieldValueType,
class ...inputFieldProperties>
1248 const Kokkos::DynRankView<inputDataValueType, inputDataProperties...> inputData,
1249 const Kokkos::DynRankView<inputFieldValueType, inputFieldProperties...> inputFields,
1250 const char transpose ) {
1252#ifdef HAVE_INTREPID2_DEBUG
1261 INTREPID2_TEST_FOR_EXCEPTION( inputData.rank() < 2 ||
1262 inputData.rank() > 4, std::invalid_argument,
1263 ">>> ERROR (ArrayTools::matmatProductDataField): inputData must have rank 2,3 or 4" );
1264 if (inputData.rank() > 2) {
1265 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(2) < 1 ||
1266 inputData.extent(2) > 3, std::invalid_argument,
1267 ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension (2) must be 1,2 or 3" );
1269 if (inputData.rank() == 4) {
1270 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(3) < 1 ||
1271 inputData.extent(3) > 3, std::invalid_argument,
1272 ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension (3) must be 1,2 or 3" );
1276 INTREPID2_TEST_FOR_EXCEPTION( inputFields.rank() < 4 ||
1277 inputFields.rank() > 5, std::invalid_argument,
1278 ">>> ERROR (ArrayTools::matmatProductDataField): inputFields must have rank 4 or 5");
1279 INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(inputFields.rank()-1) < 1 ||
1280 inputFields.extent(inputFields.rank()-1) > 3, std::invalid_argument,
1281 ">>> ERROR (ArrayTools::matmatProductDataField): inputFields dimension (rank-1) must be 1,2 or 3" );
1282 INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(inputFields.rank()-2) < 1 ||
1283 inputFields.extent(inputFields.rank()-2) > 3, std::invalid_argument,
1284 ">>> ERROR (ArrayTools::matmatProductDataField): inputFields dimension (rank-2) must be 1,2 or 3" );
1287 INTREPID2_TEST_FOR_EXCEPTION( outputFields.rank() != 5, std::invalid_argument,
1288 ">>> ERROR (ArrayTools::matmatProductDataField): outputFields must have rank 5" );
1289 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(3) < 1 ||
1290 outputFields.extent(3) > 3, std::invalid_argument,
1291 ">>> ERROR (ArrayTools::matmatProductDataField): outputFields dimension (3) must be 1,2 or 3" );
1292 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(4) < 1 ||
1293 outputFields.extent(4) > 3, std::invalid_argument,
1294 ">>> ERROR (ArrayTools::matmatProductDataField): outputFields dimension (4) must be 1,2 or 3" );
1307 if (inputData.extent(1) > 1) {
1308 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(2) != inputData.extent(1), std::invalid_argument,
1309 ">>> ERROR (ArrayTools::matmatProductDataField): outputFields dimension (2) does not match to inputData dimension (1)" );
1311 if (inputData.rank() == 2) {
1312 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(0) != inputData.extent(0), std::invalid_argument,
1313 ">>> ERROR (ArrayTools::matmatProductDataField): outputFields dimension (0) does not match to inputData dimension (0)" );
1315 if (inputData.rank() == 3) {
1316 const size_type f1[] = { 0, 3, 4 }, f2[] = { 0, 2, 2 };
1317 for (size_type i=0; i<3; ++i) {
1318 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputData.extent(f2[i]), std::invalid_argument,
1319 ">>> ERROR (ArrayTools::matmatProductDataField): outputFields dimension does not match to inputData dimension" );
1322 if (inputData.rank() == 4) {
1323 const size_type f1[] = { 0, 3, 4 }, f2[] = { 0, 2, 3 };
1324 for (size_type i=0; i<3; ++i) {
1325 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputData.extent(f2[i]), std::invalid_argument,
1326 ">>> ERROR (ArrayTools::matmatProductDataField): outputFields dimension does not match to inputData dimension" );
1333 if (inputFields.rank() == 5) {
1335 if (inputData.extent(1) > 1) {
1336 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(1) != inputFields.extent(2), std::invalid_argument,
1337 ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension (1) does not match to inputFields dimension (2)" );
1339 if (inputData.rank() == 2) {
1340 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(0) != inputFields.extent(0), std::invalid_argument,
1341 ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension (0) does not match to inputFields dimension (0)" );
1343 if (inputData.rank() == 3) {
1345 const size_type f1[] = { 0, 2, 2 }, f2[] = { 0, 3, 4 };
1346 for (size_type i=0; i<3; ++i) {
1347 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
1348 ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension does not match to inputFields dimension" );
1351 if (inputData.rank() == 4) {
1352 const size_type f1[] = { 0, 2, 3 }, f2[] = { 0, 3, 4 };
1353 for (size_type i=0; i<3; ++i) {
1354 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
1355 ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension does not match to inputFields dimension" );
1360 for (size_type i=0; i<outputFields.rank(); ++i) {
1361 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(i) != inputFields.extent(i), std::invalid_argument,
1362 ">>> ERROR (ArrayTools::matmatProductDataField): outputFields dimension does not match to inputFields dimension" );
1366 if (inputData.extent(1) > 1) {
1367 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(1) != inputFields.extent(1), std::invalid_argument,
1368 ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension (1) does not match to inputFields dimension (1)" );
1370 if (inputData.rank() == 3) {
1371 const size_type f1[] = { 2, 2 }, f2[] = { 2, 3 };
1372 for (size_type i=0; i<2; ++i) {
1373 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
1374 ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension does not match to inputFields dimension" );
1377 if (inputData.rank() == 4) {
1378 const size_type f1[] = { 2, 3 }, f2[] = { 2, 3 };
1379 for (size_type i=0; i<2; ++i) {
1380 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
1381 ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension does not match to inputFields dimension" );
1387 const size_type f1[] = { 1, 2, 3, 4 }, f2[] = { 0, 1, 2, 3 };
1388 for (size_type i=0; i<4; ++i) {
1389 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
1390 ">>> ERROR (ArrayTools::matmatProductDataField): outputFields dimension does not match to inputFields dimension" );
1401 transpose ==
't' || transpose ==
'T' );
1406 template<
typename DeviceType>
1407 template<
typename outputDataValueType,
class ...outputDataProperties,
1408 typename inputDataLeftValueType,
class ...inputDataLeftProperties,
1409 typename inputDataRightValueType,
class ...inputDataRightProperties>
1412 const Kokkos::DynRankView<inputDataLeftValueType, inputDataLeftProperties...> inputDataLeft,
1413 const Kokkos::DynRankView<inputDataRightValueType,inputDataRightProperties...> inputDataRight,
1414 const char transpose ) {
1416#ifdef HAVE_INTREPID2_DEBUG
1425 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.rank() < 2 ||
1426 inputDataLeft.rank() > 4, std::invalid_argument,
1427 ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft must have rank 2,3 or 4" );
1428 if (inputDataLeft.rank() > 2) {
1429 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(2) < 1 ||
1430 inputDataLeft.extent(2) > 3, std::invalid_argument,
1431 ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft dimension (2) must be 1,2 or 3" );
1433 if (inputDataLeft.rank() == 4) {
1434 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(3) < 1 ||
1435 inputDataLeft.extent(3) > 3, std::invalid_argument,
1436 ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft dimension (3) must be 1,2 or 3" );
1440 INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.rank() < 3 ||
1441 inputDataRight.rank() > 4, std::invalid_argument,
1442 ">>> ERROR (ArrayTools::matmatProductDataData): inputDataRight must have rank 3 or 4" );
1443 INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.extent(inputDataRight.rank()-1) < 1 ||
1444 inputDataRight.extent(inputDataRight.rank()-1) > 3, std::invalid_argument,
1445 ">>> ERROR (ArrayTools::matmatProductDataData): inputDataRight (rank-1) must be 1,2 or 3" );
1446 INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.extent(inputDataRight.rank()-2) < 1 ||
1447 inputDataRight.extent(inputDataRight.rank()-2) > 3, std::invalid_argument,
1448 ">>> ERROR (ArrayTools::matmatProductDataData): inputDataRight (rank-2) must be 1,2 or 3" );
1451 INTREPID2_TEST_FOR_EXCEPTION( outputData.rank() != 4, std::invalid_argument,
1452 ">>> ERROR (ArrayTools::matmatProductDataData): outputData must have rank 4" );
1453 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(2) < 1 ||
1454 outputData.extent(2) > 3, std::invalid_argument,
1455 ">>> ERROR (ArrayTools::matmatProductDataData): outputData (2) must be 1,2 or 3" );
1456 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(3) < 1 ||
1457 outputData.extent(3) > 3, std::invalid_argument,
1458 ">>> ERROR (ArrayTools::matmatProductDataData): outputData (3) must be 1,2 or 3" );
1471 if (inputDataLeft.extent(1) > 1) {
1472 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(1) != inputDataLeft.extent(1), std::invalid_argument,
1473 ">>> ERROR (ArrayTools::matmatProductDataData): outputData dimension (1) does not match to inputDataLeft dimension (1)" );
1475 if (inputDataLeft.rank() == 2) {
1476 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(0) != inputDataLeft.extent(0), std::invalid_argument,
1477 ">>> ERROR (ArrayTools::matmatProductDataData): outputData dimension (0) does not match to inputDataLeft dimension (0)" );
1479 if (inputDataLeft.rank() == 3) {
1480 const size_type f1[] = { 0, 2, 3 }, f2[] = { 0, 2, 2 };
1481 for (size_type i=0; i<3; ++i) {
1482 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataLeft.extent(f2[i]), std::invalid_argument,
1483 ">>> ERROR (ArrayTools::matmatProductDataData): outputData dimension does not match to inputDataLeft dimension" );
1486 if (inputDataLeft.rank() == 4) {
1487 const size_type f1[] = { 0, 2, 3 }, f2[] = { 0, 2, 3 };
1488 for (size_type i=0; i<3; ++i) {
1489 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataLeft.extent(f2[i]), std::invalid_argument,
1490 ">>> ERROR (ArrayTools::matmatProductDataData): outputData dimension does not match to inputDataLeft dimension" );
1497 if (inputDataRight.rank() == 4) {
1499 if (inputDataLeft.extent(1) > 1) {
1500 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(1) != inputDataRight.extent(1), std::invalid_argument,
1501 ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft dimension (1) does not match to inputDataRight dimension (1)" );
1503 if (inputDataLeft.rank() == 2) {
1504 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(0) != inputDataRight.extent(0), std::invalid_argument,
1505 ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft dimension (0) does not match to inputDataRight dimension (0)" );
1507 if (inputDataLeft.rank() == 3) {
1508 const size_type f1[] = { 0, 2, 2 }, f2[] = { 0, 2, 3 };
1509 for (size_type i=0; i<3; ++i) {
1510 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1511 ">>> ERROR (ArrayTools::matmatProductDataData): intputDataLeft dimension does not match to inputDataRight dimension" );
1514 if (inputDataLeft.rank() == 4) {
1515 const size_type f1[] = { 0, 2, 3 }, f2[] = { 0, 2, 3 };
1516 for (size_type i=0; i<3; ++i) {
1517 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1518 ">>> ERROR (ArrayTools::matmatProductDataData): intputDataLeft dimension does not match to inputDataRight dimension" );
1523 for (size_type i=0; i< outputData.rank(); ++i) {
1524 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(i) != inputDataRight.extent(i), std::invalid_argument,
1525 ">>> ERROR (ArrayTools::matmatProductDataData): outputData dimension does not match to inputDataRight dimension" );
1529 if (inputDataLeft.extent(1) > 1) {
1530 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(1) != inputDataRight.extent(0), std::invalid_argument,
1531 ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft dimension (1) does not match to inputDataRight dimension (0)" );
1533 if (inputDataLeft.rank() == 3) {
1534 const size_type f1[] = { 2, 2 }, f2[] = { 1, 2 };
1535 for (size_type i=0; i<2; ++i) {
1536 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1537 ">>> ERROR (ArrayTools::matmatProductDataData): intputDataLeft dimension does not match to inputDataRight dimension" );
1540 if (inputDataLeft.rank() == 4) {
1541 const size_type f1[] = { 2, 3 }, f2[] = { 1, 2 };
1542 for (size_type i=0; i<2; ++i) {
1543 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1544 ">>> ERROR (ArrayTools::matmatProductDataData): intputDataLeft dimension does not match to inputDataRight dimension" );
1549 const size_type f1[] = { 1, 2, 3 }, f2[] = { 0, 1, 2 };
1550 for (size_type i=0; i<3; ++i) {
1551 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1552 ">>> ERROR (ArrayTools::matmatProductDataData): outputData dimension does not match to inputDataRight dimension" );
1563 transpose ==
't' || transpose ==
'T' );