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