Intrepid2
Intrepid2_ArrayToolsDefDot.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Intrepid2 Package
4 //
5 // Copyright 2007 NTESS and the Intrepid2 contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
16 #ifndef __INTREPID2_ARRAYTOOLS_DEF_DOT_HPP__
17 #define __INTREPID2_ARRAYTOOLS_DEF_DOT_HPP__
18 
19 namespace Intrepid2 {
20 
21  namespace FunctorArrayTools {
25  template < typename OutputViewType, typename leftInputViewType, typename rightInputViewType >
26  struct F_dotMultiply {
27  OutputViewType _output;
28  leftInputViewType _leftInput;
29  rightInputViewType _rightInput;
30  const bool _hasField;
31  typedef typename OutputViewType::value_type value_type;
32 
33  KOKKOS_INLINE_FUNCTION
34  F_dotMultiply(OutputViewType output_,
35  leftInputViewType leftInput_,
36  rightInputViewType rightInput_,
37  const bool hasField_)
38  : _output(output_), _leftInput(leftInput_), _rightInput(rightInput_),
39  _hasField(hasField_) {}
40 
41  KOKKOS_INLINE_FUNCTION
42  void operator()(const size_type iter) const {
43  size_type cl(0), bf(0), pt(0);
44  size_type leftRank(_leftInput.rank()), rightRank(_rightInput.rank());
45 
46  if (_hasField)
47  unrollIndex( cl, bf, pt,
48  _output.extent(0),
49  _output.extent(1),
50  _output.extent(2),
51  iter );
52  else
53  unrollIndex( cl, pt,
54  _output.extent(0),
55  _output.extent(1),
56  iter);
57 
58  auto result = ( _hasField ? Kokkos::subview(_output, cl, bf, pt) :
59  Kokkos::subview(_output, cl, pt));
60 
61  const auto left = (_leftInput.extent(1) == 1) ? Kokkos::subview(_leftInput, cl, 0, Kokkos::ALL(), Kokkos::ALL()) :
62  Kokkos::subview(_leftInput, cl, pt, Kokkos::ALL(), Kokkos::ALL());
63 
64 
65  const auto right = (rightRank == leftRank + ordinal_type(_hasField)) ?
66  ( _hasField ? Kokkos::subview(_rightInput, cl, bf, pt, Kokkos::ALL(), Kokkos::ALL()) :
67  Kokkos::subview(_rightInput, cl, pt, Kokkos::ALL(), Kokkos::ALL())) :
68  ( _hasField ? Kokkos::subview(_rightInput, bf, pt, Kokkos::ALL(), Kokkos::ALL()) :
69  Kokkos::subview(_rightInput, pt, Kokkos::ALL(), Kokkos::ALL()));
70 
71  const ordinal_type iend = left.extent(0);
72  const ordinal_type jend = left.extent(1);
73 
74  value_type tmp(0);
75  for(ordinal_type i = 0; i < iend; ++i)
76  for(ordinal_type j = 0; j < jend; ++j)
77  tmp += left(i, j)*right(i, j);
78  result() = tmp;
79  }
80  };
81  } //namespace
82 
83  template<typename DeviceType>
84  template<typename outputValueType, class ...outputProperties,
85  typename leftInputValueType, class ...leftInputProperties,
86  typename rightInputValueType, class ...rightInputProperties>
87  void
89  dotMultiply( Kokkos::DynRankView<outputValueType, outputProperties...> output,
90  const Kokkos::DynRankView<leftInputValueType, leftInputProperties...> leftInput,
91  const Kokkos::DynRankView<rightInputValueType,rightInputProperties...> rightInput,
92  const bool hasField ) {
93 
94  typedef Kokkos::DynRankView<outputValueType, outputProperties...> OutputViewType;
95  typedef Kokkos::DynRankView<leftInputValueType, leftInputProperties...> leftInputViewType;
96  typedef Kokkos::DynRankView<rightInputValueType,rightInputProperties...> rightInputViewType;
98 
99  const size_type loopSize = ( hasField ? output.extent(0)*output.extent(1)*output.extent(2) :
100  output.extent(0)*output.extent(1) );
101  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
102  Kokkos::parallel_for( policy, FunctorType(output, leftInput, rightInput, hasField) );
103  }
104 
105 
106 
107  template<typename DeviceType>
108  template<typename outputFieldValueType, class ...outputFieldProperties,
109  typename inputDataValueType, class ...inputDataProperties,
110  typename inputFieldValueType, class ...inputFieldProperties>
111  void
113  dotMultiplyDataField( Kokkos::DynRankView<outputFieldValueType,outputFieldProperties...> outputFields,
114  const Kokkos::DynRankView<inputDataValueType, inputDataProperties...> inputData,
115  const Kokkos::DynRankView<inputFieldValueType, inputFieldProperties...> inputFields ) {
116 
117 #ifdef HAVE_INTREPID2_DEBUG
118  {
119  if (inputFields.rank() > inputData.rank()) {
120  INTREPID2_TEST_FOR_EXCEPTION( inputData.rank() < 2 || inputData.rank() > 4, std::invalid_argument,
121  ">>> ERROR (ArrayTools::dotMultiplyDataField): Input data container must have rank 2, 3 or 4.");
122  INTREPID2_TEST_FOR_EXCEPTION( inputFields.rank() != (inputData.rank()+1), std::invalid_argument,
123  ">>> ERROR (ArrayTools::dotMultiplyDataField): Input fields container must have rank one larger than the rank of the input data container.");
124  INTREPID2_TEST_FOR_EXCEPTION( outputFields.rank() != 3, std::invalid_argument,
125  ">>> ERROR (ArrayTools::dotMultiplyDataField): Output fields container must have rank 3.");
126  INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(0) != inputData.extent(0), std::invalid_argument,
127  ">>> ERROR (ArrayTools::dotMultiplyDataField): Zeroth dimensions (number of integration domains) of the fields and data input containers must agree!");
128  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(1) != inputFields.extent(2) &&
129  inputData.extent(1) != 1, std::invalid_argument,
130  ">>> ERROR (ArrayTools::dotMultiplyDataField): Second dimension of the fields input container and first dimension of data input container (number of integration points) must agree or first data dimension must be 1!");
131  for (size_type i=2;i<inputData.rank();++i) {
132  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(i) != inputFields.extent(i+1), std::invalid_argument,
133  ">>> ERROR (ArrayTools::dotMultiplyDataField): inputData dimension (i) does not match to the dimension (i+1) of inputFields");
134  }
135  for (size_type i=0;i<outputFields.rank();++i) {
136  INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(i) != outputFields.extent(i), std::invalid_argument,
137  ">>> ERROR (ArrayTools::dotMultiplyDataField): inputFields dimension (i) does not match to the dimension (i+1) of outputFields");
138  }
139  } else {
140  INTREPID2_TEST_FOR_EXCEPTION( inputData.rank() < 2 || inputData.rank() > 4, std::invalid_argument,
141  ">>> ERROR (ArrayTools::dotMultiplyDataField): Input data container must have rank 2, 3 or 4.");
142  INTREPID2_TEST_FOR_EXCEPTION( inputFields.rank() != inputData.rank(), std::invalid_argument,
143  ">>> ERROR (ArrayTools::dotMultiplyDataField): The rank of fields input container must equal the rank of data input container.");
144  INTREPID2_TEST_FOR_EXCEPTION( outputFields.rank() != 3, std::invalid_argument,
145  ">>> ERROR (ArrayTools::dotMultiplyDataField): Output fields container must have rank 3.");
146  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(1) != inputFields.extent(1) &&
147  inputData.extent(1) != 1, std::invalid_argument,
148  ">>> ERROR (ArrayTools::dotMultiplyDataField): First dimensions of the fields and data input containers (number of integration points) must agree or first data dimension must be 1!");
149  INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(0) != outputFields.extent(1), std::invalid_argument,
150  ">>> ERROR (ArrayTools::dotMultiplyDataField): Zeroth dimension of the fields input container and first dimension of the fields output container (number of fields) must agree!");
151  INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(1) != outputFields.extent(2), std::invalid_argument,
152  ">>> ERROR (ArrayTools::dotMultiplyDataField): First dimension of the fields input container and second dimension of the fields output container (number of integration points) must agree!");
153  INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(0) != inputData.extent(0), std::invalid_argument,
154  ">>> ERROR (ArrayTools::dotMultiplyDataField): Zeroth dimensions of the fields output and data input containers (number of integration domains) must agree!");
155  for (size_type i=2;i<inputData.rank();++i) {
156  INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(i) != inputFields.extent(i), std::invalid_argument,
157  ">>> ERROR (ArrayTools::dotMultiplyDataField): inputData dimension (i) does not match to the dimension (i) of inputFields");
158  }
159  }
160  }
161 #endif
162 
164  inputData,
165  inputFields,
166  true );
167  }
168 
169 
170 
171  template<typename DeviceType>
172  template<typename outputDataValueType, class ...outputDataProperties,
173  typename inputDataLeftValueType, class ...inputDataLeftProperties,
174  typename inputDataRightValueType, class ...inputDataRightProperties>
175  void
177  dotMultiplyDataData( Kokkos::DynRankView<outputDataValueType, outputDataProperties...> outputData,
178  const Kokkos::DynRankView<inputDataLeftValueType, inputDataLeftProperties...> inputDataLeft,
179  const Kokkos::DynRankView<inputDataRightValueType,inputDataRightProperties...> inputDataRight ) {
180 
181 #ifdef HAVE_INTREPID2_DEBUG
182  {
183  if (inputDataRight.rank() >= inputDataLeft.rank()) {
184  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.rank() < 2 || inputDataLeft.rank() > 4, std::invalid_argument,
185  ">>> ERROR (ArrayTools::dotMultiplyDataData): Left data input container must have rank 2, 3 or 4.");
186  INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.rank() != inputDataLeft.rank(), std::invalid_argument,
187  ">>> ERROR (ArrayTools::dotMultiplyDataData): The rank of the right data input container must equal the rank of the left data input container.");
188  INTREPID2_TEST_FOR_EXCEPTION( outputData.rank() != 2, std::invalid_argument,
189  ">>> ERROR (ArrayTools::dotMultiplyDataData): Data output container must have rank 2.");
190  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(1) != inputDataRight.extent(1) &&
191  inputDataLeft.extent(1) != 1, std::invalid_argument,
192  ">>> ERROR (ArrayTools::dotMultiplyDataData): First dimensions of the left and right data input containers (number of integration points) must agree or first left data dimension must be 1!");
193  for (size_type i=0;i<inputDataLeft.rank();++i) {
194  if (i != 1) {
195  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(i) != inputDataRight.extent(i), std::invalid_argument,
196  ">>> ERROR (ArrayTools::dotMultiplyDataData): inputDataLeft dimension (i) does not match to the dimension (i) of inputDataRight");
197  }
198  }
199  for (size_type i=0;i<outputData.rank();++i) {
200  INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.extent(i) != outputData.extent(i), std::invalid_argument,
201  ">>> ERROR (ArrayTools::dotMultiplyDataData): inputDataRight dimension (i) does not match to the dimension (i) of outputData");
202  }
203  } else {
204  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.rank() < 2 || inputDataLeft.rank() > 4, std::invalid_argument,
205  ">>> ERROR (ArrayTools::dotMultiplyDataData): Left data input container must have rank 2, 3 or 4.");
206  INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.rank() != (inputDataLeft.rank()-1), std::invalid_argument,
207  ">>> ERROR (ArrayTools::dotMultiplyDataData): Right data input container must have rank one less than the rank of left data input container.");
208  INTREPID2_TEST_FOR_EXCEPTION( outputData.rank() != 2, std::invalid_argument,
209  ">>> ERROR (ArrayTools::dotMultiplyDataData): Data output container must have rank 2.");
210  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(1) != inputDataRight.extent(0) &&
211  inputDataLeft.extent(1) != 1, std::invalid_argument,
212  ">>> ERROR (ArrayTools::dotMultiplyDataData): Zeroth dimension of the right data input container and first dimension of left data input container (number of integration points) must agree or first left data dimension must be 1!");
213  INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.extent(0) != outputData.extent(1), std::invalid_argument,
214  ">>> ERROR (ArrayTools::dotMultiplyDataData): Zeroth dimension of the right data input container and first dimension of output data container (number of integration points) must agree!");
215  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(0) != outputData.extent(0), std::invalid_argument,
216  ">>> ERROR (ArrayTools::dotMultiplyDataData): Zeroth dimensions of the left data input and data output containers (number of integration domains) must agree!");
217  for (size_type i=1;i<inputDataRight.rank();++i) {
218  INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(i+1) != inputDataRight.extent(i), std::invalid_argument,
219  ">>> ERROR (ArrayTools::dotMultiplyDataData): inputDataLeft dimension (i+1) does not match to the dimension (i) of inputDataRight");
220  }
221  }
222  }
223 #endif
224 
226  inputDataLeft,
227  inputDataRight,
228  false );
229  }
230 
231 } // end namespace Intrepid2
232 #endif
static void dotMultiplyDataData(Kokkos::DynRankView< outputDataValueType, outputDataProperties...> outputData, const Kokkos::DynRankView< inputDataLeftValueType, inputDataLeftProperties...> inputDataLeft, const Kokkos::DynRankView< inputDataRightValueType, inputDataRightProperties...> inputDataRight)
There are two use cases: (1) dot product of a rank-2, 3 or 4 container inputDataRight with dimensions...
Functor for dotMultiply see Intrepid2::ArrayTools for more.
Utility class that provides methods for higher-order algebraic manipulation of user-defined arrays...
static void dotMultiplyDataField(Kokkos::DynRankView< outputFieldValueType, outputFieldProperties...> outputFields, const Kokkos::DynRankView< inputDataValueType, inputDataProperties...> inputData, const Kokkos::DynRankView< inputFieldValueType, inputFieldProperties...> inputFields)
There are two use cases: (1) dot product of a rank-3, 4 or 5 container inputFields with dimensions (C...