Intrepid2
Intrepid2_ArrayToolsDefScalar.hpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Intrepid Package
5// Copyright (2007) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Kyungjoo Kim (kyukim@sandia.gov), or
38// Mauro Perego (mperego@sandia.gov)
39//
40// ************************************************************************
41// @HEADER
42
49#ifndef __INTREPID2_ARRAYTOOLS_DEF_SCALAR_HPP__
50#define __INTREPID2_ARRAYTOOLS_DEF_SCALAR_HPP__
51
52namespace Intrepid2 {
53
54 namespace FunctorArrayTools {
58 template<typename OutputViewType,
59 typename inputLeftViewType,
60 typename inputRightViewType,
61 bool equalRank,
62 bool reciprocal>
64 OutputViewType _output;
65 const inputLeftViewType _inputLeft;
66 const inputRightViewType _inputRight;
67
68 KOKKOS_INLINE_FUNCTION
69 F_scalarMultiply(OutputViewType output_,
70 inputLeftViewType inputLeft_,
71 inputRightViewType inputRight_)
72 : _output(output_),
73 _inputLeft(inputLeft_),
74 _inputRight(inputRight_) {}
75
76 // DataData
77 KOKKOS_INLINE_FUNCTION
78 void operator()(const ordinal_type cl,
79 const ordinal_type pt) const {
80 const auto val = _inputLeft(cl , pt%_inputLeft.extent(1));
81
82 //const ordinal_type cp[2] = { cl, pt };
83 //ViewAdapter<2,OutputViewType> out(cp, _output);
84 auto out = Kokkos::subview(_output, cl, pt, Kokkos::ALL(), Kokkos::ALL());
85 if (equalRank) {
86 //const ViewAdapter<2,inputRightViewType> right(cp, _inputRight);
87 const auto right = Kokkos::subview(_inputRight, cl, pt, Kokkos::ALL(), Kokkos::ALL());
88 if (reciprocal) Kernels::inv_scalar_mult_mat(out, val, right);
89 else Kernels:: scalar_mult_mat(out, val, right);
90 } else {
91 //const ordinal_type p[1] = { pt };
92 //const ViewAdapter<1,inputRightViewType> right(p, _inputRight);
93 const auto right = Kokkos::subview(_inputRight, pt, Kokkos::ALL(), Kokkos::ALL());
94 if (reciprocal) Kernels::inv_scalar_mult_mat(out, val, right);
95 else Kernels:: scalar_mult_mat(out, val, right);
96 }
97 }
98
99 // DataField
100 KOKKOS_INLINE_FUNCTION
101 void operator()(const ordinal_type cl,
102 const ordinal_type bf,
103 const ordinal_type pt) const {
104 const auto val = _inputLeft(cl , pt%_inputLeft.extent(1));
105
106 //const ordinal_type cfp[3] = { cl, bf, pt };
107 //ViewAdapter<3,OutputViewType> out(cfp, _output);
108 auto out = Kokkos::subview(_output, cl, bf, pt, Kokkos::ALL(), Kokkos::ALL());
109 if (equalRank) {
110 //const ViewAdapter<3,inputRightViewType> right(cfp, _inputRight);
111 auto right = Kokkos::subview(_inputRight, cl, bf, pt, Kokkos::ALL(), Kokkos::ALL());
112 if (reciprocal) Kernels::inv_scalar_mult_mat(out, val, right);
113 else Kernels:: scalar_mult_mat(out, val, right);
114 } else {
115 //const ordinal_type fp[2] = { bf, pt };
116 //const ViewAdapter<2,inputRightViewType> right(fp, _inputRight);
117 auto right = Kokkos::subview(_inputRight, bf, pt, Kokkos::ALL(), Kokkos::ALL());
118 if (reciprocal) Kernels::inv_scalar_mult_mat(out, val, right);
119 else Kernels:: scalar_mult_mat(out, val, right);
120 }
121 }
122 };
123 }
124
125 template<typename DeviceType>
126 template<typename outputFieldValueType, class ...outputFieldProperties,
127 typename inputDataValueType, class ...inputDataProperties,
128 typename inputFieldValueType, class ...inputFieldProperties>
129 void
131 scalarMultiplyDataField( Kokkos::DynRankView<outputFieldValueType,outputFieldProperties...> outputFields,
132 const Kokkos::DynRankView<inputDataValueType, inputDataProperties...> inputData,
133 const Kokkos::DynRankView<inputFieldValueType, inputFieldProperties...> inputFields,
134 const bool reciprocal ) {
135
136#ifdef HAVE_INTREPID2_DEBUG
137 {
138 INTREPID2_TEST_FOR_EXCEPTION( inputData.rank() != 2, std::invalid_argument,
139 ">>> ERROR (ArrayTools::scalarMultiplyDataField): Input data container must have rank 2.");
140
141 if (outputFields.rank() <= inputFields.rank()) {
142 INTREPID2_TEST_FOR_EXCEPTION( inputFields.rank() < 3 || inputFields.rank() > 5, std::invalid_argument,
143 ">>> ERROR (ArrayTools::scalarMultiplyDataField): Input fields container must have rank 3, 4, or 5.");
144 INTREPID2_TEST_FOR_EXCEPTION( outputFields.rank() != inputFields.rank(), std::invalid_argument,
145 ">>> ERROR (ArrayTools::scalarMultiplyDataField): Input and output fields containers must have the same rank.");
146 INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(0) != inputData.extent(0), std::invalid_argument,
147 ">>> ERROR (ArrayTools::scalarMultiplyDataField): Zeroth dimensions (number of integration domains) of the fields and data input containers must agree!");
148 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(1) != inputFields.extent(2) &&
149 inputData.extent(1) != 1, std::invalid_argument,
150 ">>> ERROR (ArrayTools::scalarMultiplyDataField): 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!");
151 for (size_type i=0;i<inputFields.rank();++i) {
152 INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(i) != outputFields.extent(i), std::invalid_argument,
153 ">>> ERROR (ArrayTools::scalarMultiplyDataField): inputFields dimension (i) does not match to the dimension (i) of outputFields");
154 }
155 }
156 else {
157 INTREPID2_TEST_FOR_EXCEPTION( inputFields.rank() < 2 || inputFields.rank() > 4, std::invalid_argument,
158 ">>> ERROR (ArrayTools::scalarMultiplyDataField): Input fields container must have rank 2, 3, or 4.");
159 INTREPID2_TEST_FOR_EXCEPTION( outputFields.rank() != (inputFields.rank()+1), std::invalid_argument,
160 ">>> ERROR (ArrayTools::scalarMultiplyDataField): The rank of the input fields container must be one less than the rank of the output fields container.");
161 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(1) != inputFields.extent(1) &&
162 inputData.extent(1) != 1, std::invalid_argument,
163 ">>> ERROR (ArrayTools::scalarMultiplyDataField): First dimensions of fields input container and data input container (number of integration points) must agree or first data dimension must be 1!");
164 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(0) != outputFields.extent(0), std::invalid_argument,
165 ">>> ERROR (ArrayTools::scalarMultiplyDataField): Zeroth dimensions of fields output container and data input containers (number of integration domains) must agree!");
166 for (size_type i=0;i<inputFields.rank();++i) {
167 INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(i) != outputFields.extent(i+1), std::invalid_argument,
168 ">>> ERROR (ArrayTools::scalarMultiplyDataField): inputFields dimension (i) does not match to the dimension (i+1) of outputFields");
169 }
170 }
171 }
172#endif
173
174 typedef Kokkos::DynRankView<outputFieldValueType,outputFieldProperties...> outputFieldViewType;
175 typedef Kokkos::DynRankView<inputDataValueType,inputDataProperties...> inputDataViewType;
176 typedef Kokkos::DynRankView<inputFieldValueType,inputFieldProperties...> inputFieldViewType;
177
178 using range_policy_type = Kokkos::MDRangePolicy
179 < ExecSpaceType, Kokkos::Rank<3>, Kokkos::IndexType<ordinal_type> >;
180
181 const range_policy_type policy( { 0, 0, 0 },
182 { /*C*/ outputFields.extent(0), /*F*/ outputFields.extent(1), /*P*/ outputFields.extent(2) } );
183
184 const bool equalRank = ( outputFields.rank() == inputFields.rank() );
185 if (equalRank)
186 if (reciprocal) {
188 Kokkos::parallel_for( policy, FunctorType(outputFields, inputData, inputFields) );
189 } else {
191 Kokkos::parallel_for( policy, FunctorType(outputFields, inputData, inputFields) );
192 }
193 else
194 if (reciprocal) {
196 Kokkos::parallel_for( policy, FunctorType(outputFields, inputData, inputFields) );
197 } else {
199 Kokkos::parallel_for( policy, FunctorType(outputFields, inputData, inputFields) );
200 }
201 }
202
203
204 template<typename DeviceType>
205 template<typename outputDataValueType, class ...outputDataProperties,
206 typename inputDataLeftValueType, class ...inputDataLeftProperties,
207 typename inputDataRightValueType, class ...inputDataRightProperties>
208 void
210 scalarMultiplyDataData( Kokkos::DynRankView<outputDataValueType, outputDataProperties...> outputData,
211 const Kokkos::DynRankView<inputDataLeftValueType, inputDataLeftProperties...> inputDataLeft,
212 const Kokkos::DynRankView<inputDataRightValueType,inputDataRightProperties...> inputDataRight,
213 const bool reciprocal ) {
214
215#ifdef HAVE_INTREPID2_DEBUG
216 {
217 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.rank() != 2, std::invalid_argument,
218 ">>> ERROR (ArrayTools::scalarMultiplyDataData): Left input data container must have rank 2.");
219
220 if (outputData.rank() <= inputDataRight.rank()) {
221 INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.rank() < 2 || inputDataRight.rank() > 4, std::invalid_argument,
222 ">>> ERROR (ArrayTools::scalarMultiplyDataData): Right input data container must have rank 2, 3, or 4.");
223 INTREPID2_TEST_FOR_EXCEPTION( outputData.rank() != inputDataRight.rank(), std::invalid_argument,
224 ">>> ERROR (ArrayTools::scalarMultiplyDataData): Right input and output data containers must have the same rank.");
225 INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.extent(0) != inputDataLeft.extent(0), std::invalid_argument,
226 ">>> ERROR (ArrayTools::scalarMultiplyDataData): Zeroth dimensions (number of integration domains) of the left and right data input containers must agree!");
227 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(1) != inputDataRight.extent(1) &&
228 inputDataLeft.extent(1) != 1, std::invalid_argument,
229 ">>> ERROR (ArrayTools::scalarMultiplyDataData): First dimensions of the left and right data input containers (number of integration points) must agree, or first dimension of the left data input container must be 1!");
230 for (size_type i=0;i<inputDataRight.rank();++i) {
231 INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.extent(i) != outputData.extent(i), std::invalid_argument,
232 ">>> ERROR (ArrayTools::scalarMultiplyDataData): inputDataRight dimension (i) does not match to the dimension (i) of outputData");
233 }
234 } else {
235 INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.rank() < 1 || inputDataRight.rank() > 3, std::invalid_argument,
236 ">>> ERROR (ArrayTools::scalarMultiplyDataData): Right input data container must have rank 1, 2, or 3.");
237 INTREPID2_TEST_FOR_EXCEPTION( outputData.rank() != (inputDataRight.rank()+1), std::invalid_argument,
238 ">>> ERROR (ArrayTools::scalarMultiplyDataData): The rank of the right input data container must be one less than the rank of the output data container.");
239 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(1) != inputDataRight.extent(0) &&
240 inputDataLeft.extent(1) != 1, std::invalid_argument,
241 ">>> ERROR (ArrayTools::scalarMultiplyDataData): Zeroth dimension of the right input data container and first dimension of the left data input container (number of integration points) must agree or first dimension of the left data input container must be 1!");
242 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(0) != outputData.extent(0), std::invalid_argument,
243 ">>> ERROR (ArrayTools::scalarMultiplyDataData): Zeroth dimensions of data output and left data input containers (number of integration domains) must agree!");
244 for (size_type i=0;i<inputDataRight.rank();++i) {
245 INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.extent(i) != outputData.extent(i+1), std::invalid_argument,
246 ">>> ERROR (ArrayTools::scalarMultiplyDataData): inputDataRight dimension (i) does not match to the dimension (i+1) of outputData");
247 }
248 }
249 }
250#endif
251
252 typedef Kokkos::DynRankView<outputDataValueType,outputDataProperties...> outputDataViewType;
253 typedef Kokkos::DynRankView<inputDataLeftValueType,inputDataLeftProperties...> inputDataLeftViewType;
254 typedef Kokkos::DynRankView<inputDataRightValueType,inputDataRightProperties...> inputDataRightViewType;
255
256 using range_policy_type = Kokkos::MDRangePolicy
257 < ExecSpaceType, Kokkos::Rank<2>, Kokkos::IndexType<ordinal_type> >;
258
259 const range_policy_type policy( { 0, 0 },
260 { /*C*/ outputData.extent(0), /*P*/ outputData.extent(1) } );
261
262 const bool equalRank = ( outputData.rank() == inputDataRight.rank() );
263 if (equalRank)
264 if (reciprocal) {
266 Kokkos::parallel_for( policy, FunctorType(outputData, inputDataLeft, inputDataRight) );
267 } else {
269 Kokkos::parallel_for( policy, FunctorType(outputData, inputDataLeft, inputDataRight) );
270 }
271 else
272 if (reciprocal) {
274 Kokkos::parallel_for( policy, FunctorType(outputData, inputDataLeft, inputDataRight) );
275 } else {
277 Kokkos::parallel_for( policy, FunctorType(outputData, inputDataLeft, inputDataRight) );
278 }
279 }
280
281} // end namespace Intrepid2
282
283#endif
static void scalarMultiplyDataField(Kokkos::DynRankView< outputFieldValueType, outputFieldProperties... > outputFields, const Kokkos::DynRankView< inputDataValueType, inputDataProperties... > inputData, const Kokkos::DynRankView< inputFieldValueType, inputFieldProperties... > inputFields, const bool reciprocal=false)
There are two use cases: (1) multiplies a rank-3, 4, or 5 container inputFields with dimensions (C,...
static void scalarMultiplyDataData(Kokkos::DynRankView< outputDataValueType, outputDataProperties... > outputData, const Kokkos::DynRankView< inputDataLeftValueType, inputDataLeftProperties... > inputDataLeft, const Kokkos::DynRankView< inputDataRightValueType, inputDataRightProperties... > inputDataRight, const bool reciprocal=false)
There are two use cases: (1) multiplies a rank-2, 3, or 4 container inputDataRight with dimensions (C...
Functor for scalarMultiply see Intrepid2::ArrayTools for more.