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