Intrepid2
Intrepid2_ArrayToolsDefTensor.hpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Intrepid Package
5// Copyright (2007) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Kyungjoo Kim (kyukim@sandia.gov), or
38// Mauro Perego (mperego@sandia.gov)
39//
40// ************************************************************************
41// @HEADER
42
49#ifndef __INTREPID2_ARRAYTOOLS_DEF_TENSOR_HPP__
50#define __INTREPID2_ARRAYTOOLS_DEF_TENSOR_HPP__
51
52namespace Intrepid2 {
53
54 namespace FunctorArrayTools {
58 template < typename OutputViewType, typename leftInputViewType, typename rightInputViewType >
59 struct F_crossProduct{
60 OutputViewType _output;
61 const leftInputViewType _leftInput;
62 const rightInputViewType _rightInput;
63 const bool _hasField, _isCrossProd3D;
64
65 KOKKOS_INLINE_FUNCTION
66 F_crossProduct(OutputViewType output_,
67 leftInputViewType leftInput_,
68 rightInputViewType rightInput_,
69 const bool hasField_,
70 const bool isCrossProd3D_)
71 : _output(output_), _leftInput(leftInput_), _rightInput(rightInput_),
72 _hasField(hasField_), _isCrossProd3D(isCrossProd3D_) {}
73
74 KOKKOS_INLINE_FUNCTION
75 void operator()(const size_type iter) const {
76 size_type cell, field(0), point;
77 size_type rightRank = _rightInput.rank();
78
79 if (_hasField)
80 unrollIndex( cell, field, point,
81 _output.extent(0),
82 _output.extent(1),
83 _output.extent(2),
84 iter );
85 else
86 unrollIndex( cell, point,
87 _output.extent(0),
88 _output.extent(1),
89 iter );
90
91 auto result = ( _hasField ? Kokkos::subview(_output, cell, field, point, Kokkos::ALL()) :
92 Kokkos::subview(_output, cell, point, Kokkos::ALL()));
93
94 auto left = Kokkos::subview(_leftInput, cell, point, Kokkos::ALL());
95
96 auto right = (rightRank == 2 + size_type(_hasField)) ?
97 ( _hasField ? Kokkos::subview(_rightInput, field, point, Kokkos::ALL()) :
98 Kokkos::subview(_rightInput, point, Kokkos::ALL())) :
99 ( _hasField ? Kokkos::subview(_rightInput, cell, field, point, Kokkos::ALL()) :
100 Kokkos::subview(_rightInput, cell, point, Kokkos::ALL()));
101
102 // branch prediction is possible
103 if (_isCrossProd3D) {
104 result(0) = left(1)*right(2) - left(2)*right(1);
105 result(1) = left(2)*right(0) - left(0)*right(2);
106 result(2) = left(0)*right(1) - left(1)*right(0);
107 } else {
108 result(0) = left(0)*right(1) - left(1)*right(0);
109 }
110 }
111 };
112 } //namespace
113
114 template<typename DeviceType>
115 template<typename outputValueType, class ...outputProperties,
116 typename leftInputValueType, class ...leftInputProperties,
117 typename rightInputValueType, class ...rightInputProperties>
118 void
119 ArrayTools<DeviceType>::Internal::
120 crossProduct( Kokkos::DynRankView<outputValueType, outputProperties...> output,
121 const Kokkos::DynRankView<leftInputValueType, leftInputProperties...> leftInput,
122 const Kokkos::DynRankView<rightInputValueType,rightInputProperties...> rightInput,
123 const bool hasField ) {
124
125 typedef Kokkos::DynRankView<outputValueType, outputProperties...> OutputViewType;
126 typedef const Kokkos::DynRankView<leftInputValueType, leftInputProperties...> leftInputViewType;
127 typedef const Kokkos::DynRankView<rightInputValueType,rightInputProperties...> rightInputViewType;
128 typedef FunctorArrayTools::F_crossProduct<OutputViewType, leftInputViewType, rightInputViewType> FunctorType;
129
130 const size_type loopSize = ( hasField ? output.extent(0)*output.extent(1)*output.extent(2) :
131 output.extent(0)*output.extent(1) );
132 const bool isCrossProd3D = (leftInput.extent(2) == 3);
133 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
134 Kokkos::parallel_for( policy, FunctorType(output, leftInput, rightInput, hasField, isCrossProd3D) );
135 }
136
137
138
139 template<typename DeviceType>
140 template<typename outputFieldValueType, class ...outputFieldProperties,
141 typename inputDataValueType, class ...inputDataProperties,
142 typename inputFieldValueType, class ...inputFieldProperties>
143 void
144 ArrayTools<DeviceType>::
145 crossProductDataField( Kokkos::DynRankView<outputFieldValueType,outputFieldProperties...> outputFields,
146 const Kokkos::DynRankView<inputDataValueType, inputDataProperties...> inputData,
147 const Kokkos::DynRankView<inputFieldValueType, inputFieldProperties...> inputFields ) {
148
149#ifdef HAVE_INTREPID2_DEBUG
150 {
151 /*
152 * Check array rank and spatial dimension range (if applicable)
153 * (1) inputData(C,P,D);
154 * (2) inputFields(C,F,P,D) or (F,P,D);
155 * (3) outputFields(C,F,P,D) in 3D, or (C,F,P) in 2D
156 */
157 // (1) inputData is (C, P, D) and 2 <= D <= 3 is required
158 INTREPID2_TEST_FOR_EXCEPTION( inputData.rank() != 3, std::invalid_argument,
159 ">>> ERROR (ArrayTools::crossProductDataField): inputData must have rank 3");
160 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(2) < 2 || inputData.extent(2) > 3, std::invalid_argument,
161 ">>> ERROR (ArrayTools::crossProductDataField): inputData dimension(2) must be 2 or 3");
162
163 // (2) inputFields is (C, F, P, D) or (F, P, D) and 2 <= (D=dimension(rank - 1)) <= 3 is required.
164 INTREPID2_TEST_FOR_EXCEPTION( inputFields.rank() < 3 || inputFields.rank() > 4, std::invalid_argument,
165 ">>> ERROR (ArrayTools::crossProductDataField): inputFields must have rank 3 or 4" );
166 INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(inputFields.rank()-1) < 2 ||
167 inputFields.extent(inputFields.rank()-1) > 3, std::invalid_argument,
168 ">>> ERROR (ArrayTools::crossProductDataField): inputFields dimension (rank-1) must have rank 2 or 3" );
169
170 // (3) outputFields is (C,F,P,D) in 3D and (C,F,P) in 2D => rank = inputData.extent(2) + 1
171 INTREPID2_TEST_FOR_EXCEPTION( outputFields.rank() != inputData.extent(2)+1, std::invalid_argument,
172 ">>> ERROR (ArrayTools::crossProductDataField): outputFields rank must match to inputData dimension(2)+1");
173 /*
174 * Dimension cross-checks:
175 * (1) inputData vs. inputFields
176 * (2) outputFields vs. inputData
177 * (3) outputFields vs. inputFields
178 *
179 * Cross-check (1):
180 */
181 if (inputFields.rank() == 4) {
182 // inputData(C,P,D) vs. inputFields(C,F,P,D): dimensions C, P, D must match
183 const size_type f1[] = { 0, 1, 2 }, f2[] = { 0, 2, 3 };
184 for (size_type i=0; i<3; ++i) {
185 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
186 ">>> ERROR (ArrayTools::crossProductDataField): inputData dimension does not match with inputFields");
187 }
188 } else {
189 // inputData(C,P,D) vs. inputFields(F,P,D): dimensions P, D must match
190 const size_type f1[] = { 1, 2 }, f2[] = { 1, 2 };
191 for (size_type i=0; i<2; ++i) {
192 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
193 ">>> ERROR (ArrayTools::crossProductDataField): inputData dimension does not match with inputFields");
194 }
195 }
196 /*
197 * Cross-check (2):
198 */
199 if (inputData.extent(2) == 2) {
200 // in 2D: outputFields(C,F,P) vs. inputData(C,P,D): dimensions C,P must match
201 // inputData(C,P,D) vs. inputFields(F,P,D): dimensions P, D must match
202 const size_type f1[] = { 0, 2 }, f2[] = { 0, 1 };
203 for (size_type i=0; i<2; ++i) {
204 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputData.extent(f2[i]), std::invalid_argument,
205 ">>> ERROR (ArrayTools::crossProductDataField): outputFields dimension does not match with inputData");
206 }
207 } else {
208 // in 3D: outputFields(C,F,P,D) vs. inputData(C,P,D): dimensions C,P,D must match
209 const size_type f1[] = { 0, 2, 3 }, f2[] = { 0, 1, 2 };
210 for (size_type i=0; i<3; ++i) {
211 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputData.extent(f2[i]), std::invalid_argument,
212 ">>> ERROR (ArrayTools::crossProductDataField): outputFields dimension does not match with inputData");
213 }
214 }
215 /*
216 * Cross-check (3):
217 */
218 if (inputData.extent(2) == 2) {
219 // In 2D:
220 if (inputFields.rank() == 4) {
221 // and rank-4 inputFields: outputFields(C,F,P) vs. inputFields(C,F,P,D): dimensions C,F,P must match
222 const size_type f1[] = { 0, 1, 2 }, f2[] = { 0, 1, 2 };
223 for (size_type i=0; i<3; ++i) {
224 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
225 ">>> ERROR (ArrayTools::crossProductDataField): outputFields dimension does not match with inputFields");
226 }
227 } else {
228 // and rank-3 inputFields: outputFields(C,F,P) vs. inputFields(F,P,D): dimensions F,P must match
229 const size_type f1[] = { 1, 2 }, f2[] = { 0, 1 };
230 for (size_type i=0; i<2; ++i) {
231 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
232 ">>> ERROR (ArrayTools::crossProductDataField): outputFields dimension does not match with inputFields");
233 }
234 }
235 } else {
236 // In 3D:
237 if (inputFields.rank() == 4) {
238 // and rank-4 inputFields: outputFields(C,F,P,D) vs. inputFields(C,F,P,D): all dimensions C,F,P,D must match
239 for (size_type i=0; i<outputFields.rank(); ++i) {
240 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(i) != inputFields.extent(i), std::invalid_argument,
241 ">>> ERROR (ArrayTools::crossProductDataField): outputFields dimension does not match with inputFields");
242 }
243 } else {
244 // and rank-3 inputFields: outputFields(C,F,P,D) vs. inputFields(F,P,D): dimensions F,P,D must match
245 const size_type f1[] = { 1, 2, 3 }, f2[] = { 0, 1, 2 };
246 for (size_type i=0; i<3; ++i) {
247 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
248 ">>> ERROR (ArrayTools::crossProductDataField): outputFields dimension does not match with inputFields");
249 }
250 }
251 }
252 }
253#endif
254 // body
255 ArrayTools<DeviceType>::Internal::crossProduct( outputFields,
256 inputData,
257 inputFields,
258 true );
259 }
260
261
262 template<typename DeviceType>
263 template<typename outputDataValueType, class ...outputDataProperties,
264 typename inputDataLeftValueType, class ...inputDataLeftProperties,
265 typename inputDataRightValueType, class ...inputDataRightProperties>
266 void
267 ArrayTools<DeviceType>::
268 crossProductDataData( Kokkos::DynRankView<outputDataValueType, outputDataProperties...> outputData,
269 const Kokkos::DynRankView<inputDataLeftValueType, inputDataLeftProperties...> inputDataLeft,
270 const Kokkos::DynRankView<inputDataRightValueType,inputDataRightProperties...> inputDataRight ) {
271
272#ifdef HAVE_INTREPID2_DEBUG
273 {
274 /*
275 * Check array rank and spatial dimension range (if applicable)
276 * (1) inputDataLeft(C,P,D);
277 * (2) inputDataRight(C,P,D) or (P,D);
278 * (3) outputData(C,P,D) in 3D, or (C,P) in 2D
279 */
280 // (1) inputDataLeft is (C, P, D) and 2 <= D <= 3 is required
281 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.rank() != 3, std::invalid_argument,
282 ">>> ERROR (ArrayTools::crossProductDataData): inputDataLeft must have rank 3");
283 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(2) < 2 || inputDataLeft.extent(2) > 3, std::invalid_argument,
284 ">>> ERROR (ArrayTools::crossProductDataData): inputDataLeft dimension(2) must be 2 or 3");
285
286 // (2) inputDataRight is (C, P, D) or (P, D) and 2 <= (D=dimension(rank - 1)) <= 3 is required.
287 INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.rank() < 2 || inputDataRight.rank() > 3, std::invalid_argument,
288 ">>> ERROR (ArrayTools::crossProductDataData): inputDataRight must have rank 2 or 3" );
289 INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.extent(inputDataRight.rank()-1) < 2 ||
290 inputDataRight.extent(inputDataRight.rank()-1) > 3, std::invalid_argument,
291 ">>> ERROR (ArrayTools::crossProductDataData): inputDataRight dimension (rank-1) must have rank 2 or 3" );
292
293 // (3) outputData is (C,P,D) in 3D and (C,P) in 2D => rank = inputDataLeft.extent(2)
294 INTREPID2_TEST_FOR_EXCEPTION( outputData.rank() != inputDataLeft.extent(2), std::invalid_argument,
295 ">>> ERROR (ArrayTools::crossProductDataData): outputData rank must match to inputDataLeft dimension(2)");
296 /*
297 * Dimension cross-checks:
298 * (1) inputDataLeft vs. inputDataRight
299 * (2) outputData vs. inputDataLeft
300 * (3) outputData vs. inputDataRight
301 *
302 * Cross-check (1):
303 */
304 if (inputDataRight.rank() == 3) {
305 // inputDataLeft(C,P,D) vs. inputDataRight(C,P,D): all dimensions C, P, D must match
306 for (size_type i=0; i<inputDataLeft.rank(); ++i) {
307 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(i) != inputDataRight.extent(i), std::invalid_argument,
308 ">>> ERROR (ArrayTools::crossProductDataData): inputDataLeft dimension does not match to inputDataRight");
309 }
310 }
311 // inputDataLeft(C, P,D) vs. inputDataRight(P,D): dimensions P, D must match
312 else {
313 const size_type f1[] = { 1, 2 }, f2[] = { 0, 1 };
314 for (size_type i=0; i<2; ++i) {
315 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
316 ">>> ERROR (ArrayTools::crossProductDataData): inputDataLeft dimension does not match to inputDataRight");
317 }
318 }
319 /*
320 * Cross-check (2):
321 */
322 if (inputDataLeft.extent(2) == 2) {
323 // in 2D: outputData(C,P) vs. inputDataLeft(C,P,D): dimensions C, P must match
324 const size_type f1[] = { 1, 2 }, f2[] = { 0, 1 };
325 for (size_type i=0; i<2; ++i) {
326 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != outputData.extent(f2[i]), std::invalid_argument,
327 ">>> ERROR (ArrayTools::crossProductDataData): inputDataLeft dimension does not match to outputData");
328 }
329 } else {
330 // in 3D: outputData(C,P,D) vs. inputDataLeft(C,P,D): all dimensions C, P, D must match
331 for (size_type i=0; i<inputDataLeft.rank(); ++i) {
332 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(i) != outputData.extent(i), std::invalid_argument,
333 ">>> ERROR (ArrayTools::crossProductDataData): inputDataLeft dimension does not match to outputData");
334 }
335 }
336 /*
337 * Cross-check (3):
338 */
339 if (inputDataLeft.extent(2) == 2) {
340 // In 2D:
341 if (inputDataRight.rank() == 3) {
342 // and rank-3 inputDataRight: outputData(C,P) vs. inputDataRight(C,P,D): dimensions C,P must match
343 const size_type f1[] = { 0, 1 }, f2[] = { 0, 1 };
344 for (size_type i=0; i<2; ++i) {
345 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
346 ">>> ERROR (ArrayTools::crossProductDataData): outputData dimension does not match to inputDataRight");
347 }
348 } else {
349 // and rank-2 inputDataRight: outputData(C,P) vs. inputDataRight(P,D): dimension P must match
350 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(1) != inputDataRight.extent(0), std::invalid_argument,
351 ">>> ERROR (ArrayTools::crossProductDataData): outputData dimension does not match to inputDataRight");
352 }
353 } else {
354 // In 3D:
355 if (inputDataRight.rank() == 3) {
356 // and rank-3 inputDataRight: outputData(C,P,D) vs. inputDataRight(C,P,D): all dimensions C,P,D must match
357 for (size_type i=0; i<outputData.rank(); ++i) {
358 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(i) != inputDataRight.extent(i), std::invalid_argument,
359 ">>> ERROR (ArrayTools::crossProductDataData): outputData dimension does not match to inputDataRight");
360 }
361 } else {
362 // and rank-2 inputDataRight: outputData(C,P,D) vs. inputDataRight(P,D): dimensions P, D must match
363 const size_type f1[] = { 1, 2 }, f2[] = { 0, 1 };
364 for (size_type i=0; i<2; ++i) {
365 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
366 ">>> ERROR (ArrayTools::crossProductDataData): outputData dimension does not match to inputDataRight");
367 }
368 }
369 }
370 }
371#endif
372 // body
373 ArrayTools<DeviceType>::Internal::crossProduct( outputData,
374 inputDataLeft,
375 inputDataRight,
376 false );
377 }
378
379
380 namespace FunctorArrayTools {
384 template < typename OutputViewType, typename leftInputViewType, typename rightInputViewType >
385 struct F_outerProduct {
386 OutputViewType _output;
387 const leftInputViewType _leftInput;
388 const rightInputViewType _rightInput;
389 const bool _hasField;
390
391 KOKKOS_INLINE_FUNCTION
392 F_outerProduct(OutputViewType output_,
393 leftInputViewType leftInput_,
394 rightInputViewType rightInput_,
395 const bool hasField_)
396 : _output(output_), _leftInput(leftInput_), _rightInput(rightInput_),
397 _hasField(hasField_) {}
398
399 KOKKOS_INLINE_FUNCTION
400 void operator()(const size_type iter) const {
401 size_type cell, field(0), point;
402 size_type rightRank = _rightInput.rank();
403
404 if (_hasField)
405 unrollIndex( cell, field, point,
406 _output.extent(0),
407 _output.extent(1),
408 _output.extent(2),
409 iter );
410 else
411 unrollIndex( cell, point,
412 _output.extent(0),
413 _output.extent(1),
414 iter );
415
416 auto result = ( _hasField ? Kokkos::subview(_output, cell, field, point, Kokkos::ALL(), Kokkos::ALL()) :
417 Kokkos::subview(_output, cell, point, Kokkos::ALL(), Kokkos::ALL()));
418
419 auto left = Kokkos::subview(_leftInput, cell, point, Kokkos::ALL());
420
421 auto right = (rightRank == 2 + size_type(_hasField)) ?
422 ( _hasField ? Kokkos::subview(_rightInput, field, point, Kokkos::ALL()) :
423 Kokkos::subview(_rightInput, point, Kokkos::ALL())) :
424 ( _hasField ? Kokkos::subview(_rightInput, cell, field, point, Kokkos::ALL()) :
425 Kokkos::subview(_rightInput, cell, point, Kokkos::ALL()));
426
427 const ordinal_type iend = result.extent(0);
428 const ordinal_type jend = result.extent(1);
429 for (ordinal_type i=0; i<iend; ++i)
430 for (ordinal_type j=0; j<jend; ++j)
431 result(i, j) = left(i)*right(j);
432 }
433 };
434 }
435
436 template<typename DeviceType>
437 template<typename outputValueType, class ...outputProperties,
438 typename leftInputValueType, class ...leftInputProperties,
439 typename rightInputValueType, class ...rightInputProperties>
440 void
441 ArrayTools<DeviceType>::Internal::
442 outerProduct( Kokkos::DynRankView<outputValueType, outputProperties...> output,
443 const Kokkos::DynRankView<leftInputValueType, leftInputProperties...> leftInput,
444 const Kokkos::DynRankView<rightInputValueType,rightInputProperties...> rightInput,
445 const bool hasField ) {
446
447 typedef Kokkos::DynRankView<outputValueType, outputProperties...> OutputViewType;
448 typedef const Kokkos::DynRankView<leftInputValueType, leftInputProperties...> leftInputViewType;
449 typedef const Kokkos::DynRankView<rightInputValueType,rightInputProperties...> rightInputViewType;
450 typedef FunctorArrayTools::F_outerProduct<OutputViewType, leftInputViewType, rightInputViewType> FunctorType;
451
452 const size_type loopSize = ( hasField ? output.extent(0)*output.extent(1)*output.extent(2) :
453 output.extent(0)*output.extent(1) );
454 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
455 Kokkos::parallel_for( policy, FunctorType(output, leftInput, rightInput, hasField) );
456 }
457
458
459 template<typename DeviceType>
460 template<typename outputFieldValueType, class ...outputFieldProperties,
461 typename inputDataValueType, class ...inputDataProperties,
462 typename inputFieldValueType, class ...inputFieldProperties>
463 void
464 ArrayTools<DeviceType>::
465 outerProductDataField( Kokkos::DynRankView<outputFieldValueType,outputFieldProperties...> outputFields,
466 const Kokkos::DynRankView<inputDataValueType, inputDataProperties...> inputData,
467 const Kokkos::DynRankView<inputFieldValueType, inputFieldProperties...> inputFields ) {
468
469#ifdef HAVE_INTREPID2_DEBUG
470 {
471 /*
472 * Check array rank and spatial dimension range (if applicable)
473 * (1) inputData(C,P,D);
474 * (2) inputFields(C,F,P,D) or (F,P,D);
475 * (3) outputFields(C,F,P,D,D)
476 */
477 // (1) inputData is (C, P, D) and 2 <= D <= 3 is required
478 INTREPID2_TEST_FOR_EXCEPTION( inputData.rank() != 3, std::invalid_argument,
479 ">>> ERROR (ArrayTools::outerProductDataField): inputData must have rank 3");
480 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(2) < 2 || inputData.extent(2) > 3, std::invalid_argument,
481 ">>> ERROR (ArrayTools::outerProductDataField): inputData dimension(2) must be 2 or 3");
482
483 // (2) inputFields is (C, F, P, D) or (F, P, D) and 2 <= (D=dimension(rank - 1)) <= 3 is required.
484 INTREPID2_TEST_FOR_EXCEPTION( inputFields.rank() < 3 || inputFields.rank() > 4, std::invalid_argument,
485 ">>> ERROR (ArrayTools::outerProductDataField): inputFields must have rank 3 or 4" );
486 INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(inputFields.rank()-1) < 2 ||
487 inputFields.extent(inputFields.rank()-1) < 3, std::invalid_argument,
488 ">>> ERROR (ArrayTools::outerProductDataField): inputFields dimension (rank-1) must have rank 2 or 3" );
489
490 // (3) outputFields is (C,F,P,D,D)
491 INTREPID2_TEST_FOR_EXCEPTION( outputFields.rank() != 5, std::invalid_argument,
492 ">>> ERROR (ArrayTools::outerProductDataField): outputFields must have rank 5");
493 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(3) < 2 ||
494 outputFields.extent(3) > 3, std::invalid_argument,
495 ">>> ERROR (ArrayTools::outerProductDataField): outputFields dimension(3) must be 2 or 3");
496 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(4) < 2 ||
497 outputFields.extent(4) > 3, std::invalid_argument,
498 ">>> ERROR (ArrayTools::outerProductDataField): outputFields dimension(4) must be 2 or 3");
499
500 /*
501 * Dimension cross-checks:
502 * (1) inputData vs. inputFields
503 * (2) outputFields vs. inputData
504 * (3) outputFields vs. inputFields
505 *
506 * Cross-check (2): outputFields(C,F,P,D,D) vs. inputData(C,P,D): dimensions C, P, D must match
507 */
508 {
509 const size_type f1[] = { 0, 2, 3, 4 }, f2[] = { 0, 1, 2, 2 };
510 for (size_type i=0; i<4; ++i) {
511 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputData.extent(f2[i]), std::invalid_argument,
512 ">>> ERROR (ArrayTools::outerProductDataField): outputFields dimension does not match with inputData");
513 }
514 }
515
516 /*
517 * Cross-checks (1,3):
518 */
519 if (inputFields.rank() == 4) {
520 // Cross-check (1): inputData(C,P,D) vs. inputFields(C,F,P,D): dimensions C, P, D must match
521 {
522 const size_type f1[] = { 0, 1, 2 }, f2[] = { 0, 2, 3 };
523 for (size_type i=0; i<3; ++i) {
524 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
525 ">>> ERROR (ArrayTools::outerProductDataField): inputData dimension does not match with inputFields");
526 }
527 }
528
529 // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(C,F,P,D): dimensions C, F, P, D must match
530 {
531 const size_type f1[] = { 0, 1, 2, 3, 4 }, f2[] = { 0, 1, 2, 3, 3 };
532 for (size_type i=0; i<5; ++i) {
533 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
534 ">>> ERROR (ArrayTools::outerProductDataField): outputFields dimension does not match with inputFields");
535 }
536 }
537 } else {
538 // Cross-check (1): inputData(C,P,D) vs. inputFields(F,P,D): dimensions P, D must match
539 {
540 const size_type f1[] = { 1, 2 }, f2[] = { 1, 2 };
541 for (size_type i=0; i<2; ++i) {
542 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
543 ">>> ERROR (ArrayTools::outerProductDataField): inputData dimension does not match with inputFields");
544 }
545 }
546
547 // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(F,P,D): dimensions F, P, D must match
548 {
549 const size_type f1[] = { 1, 2, 3, 4 }, f2[] = { 0, 1, 2, 2 };
550 for (size_type i=0; i<4; ++i) {
551 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
552 ">>> ERROR (ArrayTools::outerProductDataField): outputFields dimension does not match with inputFields");
553 }
554 }
555 }
556 }
557#endif
558 // body
559 ArrayTools<DeviceType>::Internal::outerProduct( outputFields,
560 inputData,
561 inputFields,
562 true );
563 }
564
565
566 template<typename DeviceType>
567 template<typename outputDataValueType, class ...outputDataProperties,
568 typename inputDataLeftValuetype, class ...inputDataLeftProperties,
569 typename inputDataRightValueType, class ...inputDataRightProperties>
570 void
571 ArrayTools<DeviceType>::
572 outerProductDataData( Kokkos::DynRankView<outputDataValueType, outputDataProperties...> outputData,
573 const Kokkos::DynRankView<inputDataLeftValuetype, inputDataLeftProperties...> inputDataLeft,
574 const Kokkos::DynRankView<inputDataRightValueType,inputDataRightProperties...> inputDataRight ) {
575
576#ifdef HAVE_INTREPID2_DEBUG
577 {
578 /*
579 * Check array rank and spatial dimension range (if applicable)
580 * (1) inputDataLeft(C,P,D);
581 * (2) inputDataRight(C,P,D) or (P,D);
582 * (3) outputData(C,P,D,D)
583 */
584 // (1) inputDataLeft is (C, P, D) and 2 <= D <= 3 is required
585 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.rank() != 3, std::invalid_argument,
586 ">>> ERROR (ArrayTools::outerProductDataField): inputDataLeft must have rank 3");
587 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(2) < 2 || inputDataLeft.extent(2) > 3, std::invalid_argument,
588 ">>> ERROR (ArrayTools::outerProductDataField): inputDataLeft dimension(2) must be 2 or 3");
589
590 // (2) inputDataRight is (C, P, D) or (P, D) and 2 <= (D=dimension(rank - 1)) <= 3 is required.
591 INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.rank() < 2 || inputDataRight.rank() > 3, std::invalid_argument,
592 ">>> ERROR (ArrayTools::outerProductDataField): inputDataRight must have rank 2 or 3" );
593 INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.extent(inputDataRight.rank()-1) < 2 ||
594 inputDataRight.extent(inputDataRight.rank()-1) < 3, std::invalid_argument,
595 ">>> ERROR (ArrayTools::outerProductDataField): inputDataRight dimension (rank-1) must have rank 2 or 3" );
596
597 // (3) outputData is (C,P,D,D)
598 INTREPID2_TEST_FOR_EXCEPTION( outputData.rank() != 4, std::invalid_argument,
599 ">>> ERROR (ArrayTools::outerProductDataField): outputData must have rank 5");
600 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(2) < 2 ||
601 outputData.extent(2) > 3, std::invalid_argument,
602 ">>> ERROR (ArrayTools::outerProductDataField): outputData dimension(3) must be 2 or 3");
603 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(3) < 2 ||
604 outputData.extent(3) > 3, std::invalid_argument,
605 ">>> ERROR (ArrayTools::outerProductDataField): outputData dimension(4) must be 2 or 3");
606
607 /*
608 * Dimension cross-checks:
609 * (1) inputDataLeft vs. inputDataRight
610 * (2) outputData vs. inputDataLeft
611 * (3) outputData vs. inputDataRight
612 *
613 * Cross-check (2): outputData(C,P,D,D) vs. inputDataLeft(C,P,D): dimensions C, P, D must match
614 */
615 {
616 const size_type f1[] = { 0, 1, 2, 3 }, f2[] = { 0, 1, 2, 2 };
617 for (size_type i=0; i<4; ++i) {
618 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataLeft.extent(f2[i]), std::invalid_argument,
619 ">>> ERROR (ArrayTools::outerProductDataField): outputData dimension does not match with inputDataLeft");
620 }
621 }
622
623 /*
624 * Cross-checks (1,3):
625 */
626 if (inputDataRight.rank() == 3) {
627 // Cross-check (1): inputDataLeft(C,P,D) vs. inputDataRight(C,P,D): all dimensions C, P, D must match
628 for (size_type i=0; i<inputDataLeft.rank(); ++i) {
629 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(i) != inputDataRight.extent(i), std::invalid_argument,
630 ">>> ERROR (ArrayTools::outerProductDataField): inputDataLeft dimension does not match with inputDataRight");
631 }
632
633 // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(C,P,D): dimensions C, P, D must match
634 {
635 const size_type f1[] = { 0, 1, 2, 3 }, f2[] = { 0, 1, 2, 2 };
636 for (size_type i=0; i<4; ++i) {
637 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
638 ">>> ERROR (ArrayTools::outerProductDataField): outputData dimension does not match with inputDataRight");
639 }
640 }
641 } else {
642 // Cross-check (1): inputDataLeft(C,P,D) vs. inputDataRight(P,D): dimensions P, D must match
643 {
644 const size_type f1[] = { 1, 2 }, f2[] = { 0, 1 };
645 for (size_type i=0; i<2; ++i) {
646 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
647 ">>> ERROR (ArrayTools::outerProductDataField): inputDataLeft dimension does not match with inputDataRight");
648 }
649 }
650
651 // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(P,D): dimensions P, D must match
652 {
653 const size_type f1[] = { 1, 2, 3 }, f2[] = { 0, 1, 1 };
654 for (size_type i=0; i<3; ++i) {
655 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
656 ">>> ERROR (ArrayTools::outerProductDataField): outputData dimension does not match with inputDataRight");
657 }
658 }
659 }
660 }
661#endif
662 // body
663 ArrayTools<DeviceType>::Internal::outerProduct( outputData,
664 inputDataLeft,
665 inputDataRight,
666 false );
667 }
668
669
670 namespace FunctorArrayTools {
674 template < typename OutputViewType,
675 typename leftInputViewType,
676 typename rightInputViewType>
677 struct F_matvecProduct {
678 OutputViewType _output;
679 const leftInputViewType _leftInput;
680 const rightInputViewType _rightInput;
681
682 const bool _isTranspose;
683
684 KOKKOS_INLINE_FUNCTION
685 F_matvecProduct(OutputViewType output_,
686 leftInputViewType leftInput_,
687 rightInputViewType rightInput_,
688 const bool isTranspose_)
689 : _output(output_), _leftInput(leftInput_), _rightInput(rightInput_), _isTranspose(isTranspose_) {}
690
691 template<typename resultViewType,
692 typename leftViewType,
693 typename rightViewType>
694 KOKKOS_FORCEINLINE_FUNCTION
695 static void
696 apply_matvec_product( resultViewType &result,
697 const leftViewType &left,
698 const rightViewType &right,
699 const bool isTranspose) {
700 const ordinal_type iend = result.extent(0);
701 const ordinal_type jend = right.extent(0);
702
703 typedef typename resultViewType::value_type value_type;
704
705 switch (left.rank()) {
706 case 2:
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 += left(j, i)*right(j);
712 result(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 += left(i, j)*right(j);
719 result(i) = tmp;
720 }
721 }
722 break;
723 case 1: { //matrix is diagonal
724 for (ordinal_type i=0;i<iend;++i)
725 result(i) = left(i)*right(i);
726 break;
727 }
728 case 0: { //matrix is a scaled identity
729 const value_type val = left();
730 for (ordinal_type i=0;i<iend;++i) {
731 result(i) = val*right(i);
732 }
733 break;
734 }
735 }
736 }
737
738 KOKKOS_INLINE_FUNCTION
739 void operator()(const ordinal_type cl,
740 const ordinal_type pt) const {
741 const auto rightRank = _rightInput.rank();
742 const auto leftRank = _leftInput.rank();
743
744 auto result = Kokkos::subview(_output, cl, pt, Kokkos::ALL());
745
746 const auto lpt = (_leftInput.extent(1) == 1 ? size_type(0) : pt);
747 const auto left = ( leftRank == 4 ? Kokkos::subview(_leftInput, cl, lpt, Kokkos::ALL(), Kokkos::ALL()) :
748 leftRank == 3 ? Kokkos::subview(_leftInput, cl, lpt, Kokkos::ALL()) :
749 Kokkos::subview(_leftInput, cl, lpt));
750
751 const auto right = ( rightRank == 2 ? Kokkos::subview(_rightInput, pt, Kokkos::ALL()) :
752 Kokkos::subview(_rightInput, cl, pt, Kokkos::ALL()) );
753 apply_matvec_product( result, left, right, _isTranspose );
754 }
755
756 KOKKOS_INLINE_FUNCTION
757 void operator()(const ordinal_type cl,
758 const ordinal_type bf,
759 const ordinal_type pt) const {
760 const auto rightRank = _rightInput.rank();
761 const auto leftRank = _leftInput.rank();
762
763 auto result = Kokkos::subview(_output, cl, bf, pt, Kokkos::ALL());
764
765 const auto lpt = (_leftInput.extent(1) == 1 ? size_type(0) : pt);
766 const auto left = ( leftRank == 4 ? Kokkos::subview(_leftInput, cl, lpt, Kokkos::ALL(), Kokkos::ALL()) :
767 leftRank == 3 ? Kokkos::subview(_leftInput, cl, lpt, Kokkos::ALL()) :
768 Kokkos::subview(_leftInput, cl, lpt));
769
770 const auto right = ( rightRank == 3 ? Kokkos::subview(_rightInput, bf, pt, Kokkos::ALL()) :
771 Kokkos::subview(_rightInput, cl, bf, pt, Kokkos::ALL()));
772
773 apply_matvec_product( result, left, right, _isTranspose );
774 }
775 };
776 } //namespace
777
778 template<typename DeviceType>
779 template<typename outputValueType, class ...outputProperties,
780 typename leftInputValueType, class ...leftInputProperties,
781 typename rightInputValueType, class ...rightInputProperties>
782 void
783 ArrayTools<DeviceType>::Internal::
784 matvecProduct( Kokkos::DynRankView<outputValueType, outputProperties...> output,
785 const Kokkos::DynRankView<leftInputValueType, leftInputProperties...> leftInput,
786 const Kokkos::DynRankView<rightInputValueType,rightInputProperties...> rightInput,
787 const bool hasField,
788 const bool isTranspose ) {
789
790 typedef Kokkos::DynRankView<outputValueType, outputProperties...> OutputViewType;
791 typedef const Kokkos::DynRankView<leftInputValueType, leftInputProperties...> leftInputViewType;
792 typedef const Kokkos::DynRankView<rightInputValueType,rightInputProperties...> rightInputViewType;
793 typedef FunctorArrayTools::F_matvecProduct<OutputViewType, leftInputViewType, rightInputViewType> FunctorType;
794
795 if (hasField) {
796 using range_policy_type = Kokkos::Experimental::MDRangePolicy
797 < ExecSpaceType, Kokkos::Experimental::Rank<3>, Kokkos::IndexType<ordinal_type> >;
798 range_policy_type policy( { 0, 0, 0 },
799 { output.extent(0), output.extent(1), output.extent(2) } );
800 Kokkos::parallel_for( policy, FunctorType(output, leftInput, rightInput, isTranspose) );
801 } else {
802 using range_policy_type = Kokkos::Experimental::MDRangePolicy
803 < ExecSpaceType, Kokkos::Experimental::Rank<2>, Kokkos::IndexType<ordinal_type> >;
804 range_policy_type policy( { 0, 0 },
805 { output.extent(0), output.extent(1) } );
806 Kokkos::parallel_for( policy, FunctorType(output, leftInput, rightInput, isTranspose) );
807 }
808 }
809
810 template<typename DeviceType>
811 template<typename outputFieldValueType, class ...outputFieldProperties,
812 typename inputDataValueType, class ...inputDataProperties,
813 typename inputFieldValueType, class ...inputFieldProperties>
814 void
815 ArrayTools<DeviceType>::matvecProductDataField( Kokkos::DynRankView<outputFieldValueType,outputFieldProperties...> outputFields,
816 const Kokkos::DynRankView<inputDataValueType, inputDataProperties...> inputData,
817 const Kokkos::DynRankView<inputFieldValueType, inputFieldProperties...> inputFields,
818 const char transpose ) {
819
820#ifdef HAVE_INTREPID2_DEBUG
821 {
822 /*
823 * Check array rank and spatial dimension range (if applicable)
824 * (1) inputData(C,P), (C,P,D) or (C,P,D,D); P=1 is admissible to allow multiply by const. data
825 * (2) inputFields(C,F,P,D) or (F,P,D);
826 * (3) outputFields(C,F,P,D)
827 */
828 // (1) inputData is (C,P), (C, P, D) or (C, P, D, D) and 1 <= D <= 3 is required
829 INTREPID2_TEST_FOR_EXCEPTION( inputData.rank() < 2 || inputData.rank() > 4, std::invalid_argument,
830 ">>> ERROR (ArrayTools::matvecProductDataField): inputData must have rank 2,3 or 4" );
831 if (inputData.rank() > 2) {
832 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(2) < 1 ||
833 inputData.extent(2) > 3, std::invalid_argument,
834 ">>> ERROR (ArrayTools::matvecProductDataField): inputData dimension(2) must be 1,2 or 3");
835 }
836 if (inputData.rank() == 4) {
837 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(3) < 1 ||
838 inputData.extent(3) > 3, std::invalid_argument,
839 ">>> ERROR (ArrayTools::matvecProductDataField): inputData dimension(3) must be 1,2 or 3");
840 }
841
842 // (2) inputFields is (C, F, P, D) or (F, P, D) and 1 <= (D=dimension(rank - 1)) <= 3 is required.
843 INTREPID2_TEST_FOR_EXCEPTION( inputFields.rank() < 3 ||
844 inputFields.rank() > 4, std::invalid_argument,
845 ">>> ERROR (ArrayTools::matvecProductDataField): inputFields must have rank 3 or 4" );
846 INTREPID2_TEST_FOR_EXCEPTION( (inputFields.rank()-1) < 1 ||
847 (inputFields.rank()-1) > 3, std::invalid_argument,
848 ">>> ERROR (ArrayTools::matvecProductDataField): inputFields dimension(rank-1) must be 1,2, or 3" );
849
850 // (3) outputFields is (C,F,P,D)
851 INTREPID2_TEST_FOR_EXCEPTION( outputFields.rank() != 4, std::invalid_argument,
852 ">>> ERROR (ArrayTools::matvecProductDataField): outputFields must have rank 4" );
853 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(3) < 1 ||
854 outputFields.extent(3) > 3, std::invalid_argument,
855 ">>> ERROR (ArrayTools::matvecProductDataField): outputFields dimension(3) must be 1,2 or 3" );
856
857 /*
858 * Dimension cross-checks:
859 * (1) inputData vs. inputFields
860 * (2) outputFields vs. inputData
861 * (3) outputFields vs. inputFields
862 *
863 * Cross-check (2): outputFields(C,F,P,D) vs. inputData(C,P), (C,P,D) or (C,P,D,D):
864 * dimensions C, and D must match in all cases, dimension P must match only when non-constant
865 * data is specified (P>1). Do not check P dimensions with constant data, i.e., when P=1 in
866 * inputData(C,1,...)
867 */
868 if (inputData.extent(1) > 1) { // check P dimension if P>1 in inputData
869 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(2) != inputData.extent(1), std::invalid_argument,
870 ">>> ERROR (ArrayTools::matvecProductDataField): outputFields dimension(2) must match to inputData dimension(1)" );
871 }
872 if (inputData.rank() == 2) { // inputData(C,P) -> C match
873 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(0) != inputData.extent(0), std::invalid_argument,
874 ">>> ERROR (ArrayTools::matvecProductDataField): outputFields dimension(0) must match to inputData dimension(0)" );
875 }
876 if (inputData.rank() == 3) { // inputData(C,P,D) -> C, D match
877 const size_type f1[] = { 0, 3 }, f2[] = { 0, 2 };
878 for (size_type i=0; i<2; ++i) {
879 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputData.extent(f2[i]), std::invalid_argument,
880 ">>> ERROR (ArrayTools::matvecProductDataField): outputFields dimension does not match to inputData dimension" );
881 }
882 }
883 if (inputData.rank() == 4) { // inputData(C,P,D,D) -> C, D, D match
884 const size_type f1[] = { 0, 3, 3 }, f2[] = { 0, 2, 3 };
885 for (size_type i=0; i<3; ++i) {
886 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputData.extent(f2[i]), std::invalid_argument,
887 ">>> ERROR (ArrayTools::matvecProductDataField): outputFields dimension does not match to inputData dimension" );
888 }
889 }
890
891 /*
892 * Cross-checks (1,3):
893 */
894 if (inputFields.rank() == 4) {
895 // 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
896 if (inputData.extent(1) > 1) { // check P dimension if P>1 in inputData
897 INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(2) != inputData.extent(1), std::invalid_argument,
898 ">>> ERROR (ArrayTools::matvecProductDataField): inputFields dimension (2) does not match to inputData dimension (1)" );
899 }
900 if (inputData.rank() == 2) { // inputData(C,P) -> C match
901 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(0) != inputFields.extent(0), std::invalid_argument,
902 ">>> ERROR (ArrayTools::matvecProductDataField): inputData dimension (0) does not match to inputFields dimension (0)" );
903 }
904 if (inputData.rank() == 3) { // inputData(C,P,D) -> C, D match
905 const size_type f1[] = { 0, 2 }, f2[] = { 0, 3 };
906 for (size_type i=0; i<2; ++i) {
907 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
908 ">>> ERROR (ArrayTools::matvecProductDataField): inputData dimension does not match to inputFields dimension" );
909 }
910 }
911 if (inputData.rank() == 4) { // inputData(C,P,D,D) -> C, D, D match
912 const size_type f1[] = { 0, 2, 3 }, f2[] = { 0, 3, 3 };
913 for (size_type i=0; i<3; ++i) {
914 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
915 ">>> ERROR (ArrayTools::matvecProductDataField): inputData dimension does not match to inputFields dimension" );
916 }
917 }
918
919 // Cross-check (3): outputFields(C,F,P,D) vs. inputFields(C,F,P,D): all dimensions C, F, P, D must match
920 for (size_type i=0; i<outputFields.rank(); ++i) {
921 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(i) != inputFields.extent(i), std::invalid_argument,
922 ">>> ERROR (ArrayTools::matvecProductDataField): outputFields dimension does not match to inputFields dimension" );
923 }
924 } else {
925 // Cross-check (1): inputData(C,P), (C,P,D) or (C,P,D,D) vs. inputFields(F,P,D): dimensions P, D must match
926 if (inputData.extent(1) > 1) { // check P if P>1 in inputData
927 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(1) != inputFields.extent(1), std::invalid_argument,
928 ">>> ERROR (ArrayTools::matvecProductDataField): inputData dimension(1) does not match to inputFields dimension (1)" );
929 }
930 if (inputData.rank() == 3) {
931 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(2) != inputFields.extent(2), std::invalid_argument,
932 ">>> ERROR (ArrayTools::matvecProductDataField): inputData dimension(2) does not match to inputFields dimension (2)" );
933 }
934 if (inputData.rank() == 4) {
935 const size_type f1[] = { 2, 3 }, f2[] = { 2, 2 };
936 for (size_type i=0; i<2; ++i) {
937 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
938 ">>> ERROR (ArrayTools::matvecProductDataField): inputData dimension does not match to inputFields dimension" );
939 }
940 }
941
942 // Cross-check (3): outputFields(C,F,P,D) vs. inputFields(F,P,D): dimensions F, P, D must match
943 {
944 const size_type f1[] = { 1, 2, 3 }, f2[] = { 0, 1, 2 };
945 for (size_type i=0; i<3; ++i) {
946 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
947 ">>> ERROR (ArrayTools::matvecProductDataField): outputFields dimension does not match to inputFields dimension" );
948 }
949 }
950 }
951 }
952#endif
953 // body
954 ArrayTools<DeviceType>::Internal::matvecProduct( outputFields,
955 inputData,
956 inputFields,
957 true,
958 transpose == 't' || transpose == 'T' );
959 }
960
961
962
963 template<typename DeviceType>
964 template<typename outputDataValueType, class ...outputDataProperties,
965 typename inputDataLeftValueType, class ...inputDataLeftProperties,
966 typename inputDataRightValueType, class ...inputDataRightProperties>
967 void
968 ArrayTools<DeviceType>::
969 matvecProductDataData( Kokkos::DynRankView<outputDataValueType, outputDataProperties...> outputData,
970 const Kokkos::DynRankView<inputDataLeftValueType, inputDataLeftProperties...> inputDataLeft,
971 const Kokkos::DynRankView<inputDataRightValueType,inputDataRightProperties...> inputDataRight,
972 const char transpose ) {
973
974#ifdef HAVE_INTREPID2_DEBUG
975 {
976 /*
977 * Check array rank and spatial dimension range (if applicable)
978 * (1) inputDataLeft(C,P), (C,P,D) or (C,P,D,D); P=1 is admissible to allow multiply by const. left data
979 * (2) inputDataRight(C,P,D) or (P,D);
980 * (3) outputData(C,P,D)
981 */
982 // (1) inputDataLeft is (C,P), (C,P,D) or (C,P,D,D) and 1 <= D <= 3 is required
983 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.rank() < 2 ||
984 inputDataLeft.rank() > 4, std::invalid_argument,
985 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft must have rank 2,3 or 4" );
986
987 if (inputDataLeft.rank() > 2) {
988 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(2) < 1 ||
989 inputDataLeft.extent(2) > 3, std::invalid_argument,
990 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension(2) must be 1, 2 or 3");
991 }
992 if (inputDataLeft.rank() == 4) {
993 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(3) < 1 ||
994 inputDataLeft.extent(3) > 3, std::invalid_argument,
995 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension(3) must be 1, 2 or 3");
996 }
997
998 // (2) inputDataRight is (C, P, D) or (P, D) and 1 <= (D=dimension(rank - 1)) <= 3 is required.
999 INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.rank() < 2 ||
1000 inputDataRight.rank() > 3, std::invalid_argument,
1001 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataRight must have rank 2 or 3" );
1002 INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.extent(inputDataRight.rank()-1) < 1 ||
1003 inputDataRight.extent(inputDataRight.rank()-1) > 3, std::invalid_argument,
1004 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataRight dimension (rank-1) must be 1,2 or 3" );
1005
1006 // (3) outputData is (C,P,D)
1007 INTREPID2_TEST_FOR_EXCEPTION( outputData.rank() != 3, std::invalid_argument,
1008 ">>> ERROR (ArrayTools::matvecProductDataData): outputData must have rank 3" );
1009 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(2) < 1 ||
1010 outputData.extent(2) > 3, std::invalid_argument,
1011 ">>> ERROR (ArrayTools::matvecProductDataData): outputData dimension(2) must be 1, 2 or 3");
1012
1013 /*
1014 * Dimension cross-checks:
1015 * (1) inputDataLeft vs. inputDataRight
1016 * (2) outputData vs. inputDataLeft
1017 * (3) outputData vs. inputDataRight
1018 *
1019 * Cross-check (2): outputData(C,P,D) vs. inputDataLeft(C,P), (C,P,D) or (C,P,D,D):
1020 * dimensions C, and D must match in all cases, dimension P must match only when non-constant
1021 * data is specified (P>1). Do not check P dimensions with constant left data, i.e., when P=1 in
1022 * inputDataLeft(C,1,...)
1023 */
1024 if (inputDataLeft.extent(1) > 1) { // check P dimension if P>1 in inputDataLeft
1025 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(1) != inputDataLeft.extent(1), std::invalid_argument,
1026 ">>> ERROR (ArrayTools::matvecProductDataData): outputData dimension(1) muat match inputDataLeft dimension(1)" );
1027 }
1028 if (inputDataLeft.rank() == 2) { // inputDataLeft(C,P): check C
1029 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(0) != inputDataLeft.extent(0), std::invalid_argument,
1030 ">>> ERROR (ArrayTools::matvecProductDataData): outputData dimension(0) muat match inputDataLeft dimension(0)" );
1031 }
1032 if (inputDataLeft.rank() == 3) { // inputDataLeft(C,P,D): check C and D
1033 const size_type f1[] = { 0, 2 }, f2[] = { 0, 2 };
1034 for (size_type i=0; i<2; ++i) {
1035 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataLeft.extent(f2[i]), std::invalid_argument,
1036 ">>> ERROR (ArrayTools::matvecProductDataData): outputData dimension muat match inputDataLeft dimension" );
1037 }
1038 }
1039 if (inputDataLeft.rank() == 4) { // inputDataLeft(C,P,D,D): check C and D
1040 const size_type f1[] = { 0, 2, 2 }, f2[] = { 0, 2, 3 };
1041 for (size_type i=0; i<3; ++i) {
1042 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataLeft.extent(f2[i]), std::invalid_argument,
1043 ">>> ERROR (ArrayTools::matvecProductDataData): outputData dimension muat match inputDataLeft dimension" );
1044 }
1045 }
1046
1047 /*
1048 * Cross-checks (1,3):
1049 */
1050 if (inputDataRight.rank() == 3) {
1051 // 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
1052 if (inputDataLeft.extent(1) > 1) { // check P dimension if P>1 in inputDataLeft
1053 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(1) != inputDataRight.extent(1), std::invalid_argument,
1054 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension (1) muat match to inputDataRight dimension (1)" );
1055 }
1056 if (inputDataLeft.rank() == 2) { // inputDataLeft(C,P): check C
1057 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(0) != inputDataRight.extent(0), std::invalid_argument,
1058 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension (0) muat match to inputDataRight dimension (0)" );
1059 }
1060 if (inputDataLeft.rank() == 3) { // inputDataLeft(C,P,D): check C and D
1061 const size_type f1[] = { 0, 2 }, f2[] = { 0, 2 };
1062 for (size_type i=0; i<2; ++i) {
1063 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1064 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension muat match to inputDataRight dimension" );
1065 }
1066 }
1067 if (inputDataLeft.rank() == 4) { // inputDataLeft(C,P,D,D): check C and D
1068 const size_type f1[] = { 0, 2, 3 }, f2[] = { 0, 2, 2 };
1069 for (size_type i=0; i<3; ++i) {
1070 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1071 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension muat match to inputDataRight dimension" );
1072 }
1073 }
1074
1075 // Cross-check (3): outputData(C,P,D) vs. inputDataRight(C,P,D): all dimensions C, P, D must match
1076 for (size_type i=0; i<outputData.rank(); ++i) {
1077 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(i) != inputDataRight.extent(i), std::invalid_argument,
1078 ">>> ERROR (ArrayTools::matvecProductDataData): outputData dimension muat match to inputDataRight dimension" );
1079 }
1080 } else {
1081 // Cross-check (1): inputDataLeft(C,P), (C,P,D), or (C,P,D,D) vs. inputDataRight(P,D): dimensions P, D must match
1082 if (inputDataLeft.extent(1) > 1) { // check P if P>1 in inputData
1083 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(1) != inputDataRight.extent(0), std::invalid_argument,
1084 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension (1) does mot match to inputDataright dimension (0)" );
1085 }
1086 if (inputDataLeft.rank() == 3) { // inputDataLeft(C,P,D): check D
1087 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(2) != inputDataRight.extent(1), std::invalid_argument,
1088 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension (2) does mot match to inputDataright dimension (1)" );
1089 }
1090 if (inputDataLeft.rank() == 4) { // inputDataLeft(C,P,D,D): check D
1091 const size_type f1[] = { 2, 3 }, f2[] = { 1, 1 };
1092 for (size_type i=0; i<2; ++i) {
1093 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1094 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft dimension muat match to inputDataRight dimension" );
1095 }
1096 }
1097
1098 // Cross-check (3): outputData(C,P,D) vs. inputDataRight(P,D): dimensions P, D must match
1099 {
1100 const size_type f1[] = { 1, 2 }, f2[] = { 0, 1 };
1101 for (size_type i=0; i<2; ++i) {
1102 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1103 ">>> ERROR (ArrayTools::matvecProductDataData): outputData dimension muat match to inputDataRight dimension" );
1104 }
1105 }
1106 }
1107 }
1108#endif
1109 // body
1110 ArrayTools<DeviceType>::Internal::matvecProduct( outputData,
1111 inputDataLeft,
1112 inputDataRight,
1113 false,
1114 transpose == 't' || transpose == 'T' );
1115 }
1116
1117
1118 namespace FunctorArrayTools {
1122 template < typename OutputViewType, typename leftInputViewType, typename rightInputViewType >
1123 struct F_matmatProduct{
1124 OutputViewType _output;
1125 leftInputViewType _leftInput;
1126 rightInputViewType _rightInput;
1127 const bool _hasField, _isTranspose;
1128 typedef typename leftInputViewType::value_type value_type;
1129
1130 KOKKOS_INLINE_FUNCTION
1131 F_matmatProduct(OutputViewType output_,
1132 leftInputViewType leftInput_,
1133 rightInputViewType rightInput_,
1134 const bool hasField_,
1135 const bool isTranspose_)
1136 : _output(output_), _leftInput(leftInput_), _rightInput(rightInput_),
1137 _hasField(hasField_), _isTranspose(isTranspose_) {}
1138
1139 KOKKOS_INLINE_FUNCTION
1140 void operator()(const size_type iter) const {
1141 size_type cell(0), field(0), point(0);
1142 size_type leftRank = _leftInput.rank();
1143 size_type rightRank = _rightInput.rank();
1144
1145 if (_hasField)
1146 unrollIndex( cell, field, point,
1147 _output.extent(0),
1148 _output.extent(1),
1149 _output.extent(2),
1150 iter );
1151 else
1152 unrollIndex( cell, point,
1153 _output.extent(0),
1154 _output.extent(1),
1155 iter );
1156
1157 auto result = ( _hasField ? Kokkos::subview(_output, cell, field, point, Kokkos::ALL(), Kokkos::ALL()) :
1158 Kokkos::subview(_output, cell, point, Kokkos::ALL(), Kokkos::ALL()));
1159
1160 const auto lpoint = (_leftInput.extent(1) == 1 ? 0 : point);
1161 auto left = ( leftRank == 4 ? Kokkos::subview(_leftInput, cell, lpoint, Kokkos::ALL(), Kokkos::ALL()) :
1162 leftRank == 3 ? Kokkos::subview(_leftInput, cell, lpoint, Kokkos::ALL()) :
1163 Kokkos::subview(_leftInput, cell, lpoint) );
1164
1165 //temporary
1166 const bool hasCell = (_hasField ? rightRank == 5 : rightRank == 4);
1167 auto right = ( _hasField ? ( hasCell ? Kokkos::subview(_rightInput, cell, field, point, Kokkos::ALL(), Kokkos::ALL()) :
1168 Kokkos::subview(_rightInput, field, point, Kokkos::ALL(), Kokkos::ALL())):
1169 ( hasCell ? Kokkos::subview(_rightInput, cell, point, Kokkos::ALL(), Kokkos::ALL()) :
1170 Kokkos::subview(_rightInput, point, Kokkos::ALL(), Kokkos::ALL())));
1171
1172 const ordinal_type iend = result.extent(0);
1173 const ordinal_type jend = result.extent(1);
1174
1175 switch (leftRank) {
1176 case 4: {
1177 if (_isTranspose) {
1178 const size_type kend = right.extent(0);
1179 for (ordinal_type i=0; i<iend; ++i)
1180 for (ordinal_type j=0; j<jend; ++j) {
1181 result(i, j) = value_type(0);
1182 for (size_type k=0; k<kend; ++k)
1183 result(i, j) += left(k, i) * right(k, j);
1184 }
1185 } else {
1186 const size_type kend = right.extent(0);
1187 for (ordinal_type i=0; i<iend; ++i)
1188 for (ordinal_type j=0; j<jend; ++j) {
1189 result(i, j) = value_type(0);
1190 for (size_type k=0; k<kend; ++k)
1191 result(i, j) += left(i, k) * right(k, j);
1192 }
1193 }
1194 break;
1195 }
1196 case 3: { //matrix is diagonal
1197 for (ordinal_type i=0; i<iend; ++i)
1198 for (ordinal_type j=0; j<jend; ++j)
1199 result(i, j) = left(i) * right(i, j);
1200 break;
1201 }
1202 case 2: { //matrix is a scaled identity
1203 for (ordinal_type i=0; i<iend; ++i)
1204 for (ordinal_type j=0; j<jend; ++j)
1205 result(i, j) = left() * right(i, j);
1206 break;
1207 }
1208 }
1209 }
1210 };
1211 } //namespace
1212
1213 template<typename DeviceType>
1214 template<typename outputValueType, class ...outputProperties,
1215 typename leftInputValueType, class ...leftInputProperties,
1216 typename rightInputValueType, class ...rightInputProperties>
1217 void
1218 ArrayTools<DeviceType>::Internal::
1219 matmatProduct( Kokkos::DynRankView<outputValueType, outputProperties...> output,
1220 const Kokkos::DynRankView<leftInputValueType, leftInputProperties...> leftInput,
1221 const Kokkos::DynRankView<rightInputValueType,rightInputProperties...> rightInput,
1222 const bool hasField,
1223 const bool isTranspose ) {
1224
1225 typedef Kokkos::DynRankView<outputValueType, outputProperties...> OutputViewType;
1226 typedef const Kokkos::DynRankView<leftInputValueType, leftInputProperties...> leftInputViewType;
1227 typedef const Kokkos::DynRankView<rightInputValueType,rightInputProperties...> rightInputViewType;
1228 typedef FunctorArrayTools::F_matmatProduct<OutputViewType, leftInputViewType, rightInputViewType> FunctorType;
1229
1230 const size_type loopSize = ( hasField ? output.extent(0)*output.extent(1)*output.extent(2) :
1231 output.extent(0)*output.extent(1) );
1232 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
1233 Kokkos::parallel_for( policy, FunctorType(output, leftInput, rightInput, hasField, isTranspose) );
1234 }
1235
1236
1237
1238
1239 template<typename DeviceType>
1240 template<typename outputFieldValueType, class ...outputFieldProperties,
1241 typename inputDataValueType, class ...inputDataProperties,
1242 typename inputFieldValueType, class ...inputFieldProperties>
1243 void
1244 ArrayTools<DeviceType>::
1245 matmatProductDataField( Kokkos::DynRankView<outputFieldValueType,outputFieldProperties...> outputFields,
1246 const Kokkos::DynRankView<inputDataValueType, inputDataProperties...> inputData,
1247 const Kokkos::DynRankView<inputFieldValueType, inputFieldProperties...> inputFields,
1248 const char transpose ) {
1249
1250#ifdef HAVE_INTREPID2_DEBUG
1251 {
1252 /*
1253 * Check array rank and spatial dimension range (if applicable)
1254 * (1) inputData(C,P), (C,P,D) or (C,P,D,D); P=1 is admissible to allow multiply by const. data
1255 * (2) inputFields(C,F,P,D,D) or (F,P,D,D);
1256 * (3) outputFields(C,F,P,D,D)
1257 */
1258 // (1) inputData is (C,P), (C, P, D) or (C, P, D, D) and 1 <= D <= 3 is required
1259 INTREPID2_TEST_FOR_EXCEPTION( inputData.rank() < 2 ||
1260 inputData.rank() > 4, std::invalid_argument,
1261 ">>> ERROR (ArrayTools::matmatProductDataField): inputData must have rank 2,3 or 4" );
1262 if (inputData.rank() > 2) {
1263 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(2) < 1 ||
1264 inputData.extent(2) > 3, std::invalid_argument,
1265 ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension (2) must be 1,2 or 3" );
1266 }
1267 if (inputData.rank() == 4) {
1268 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(3) < 1 ||
1269 inputData.extent(3) > 3, std::invalid_argument,
1270 ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension (3) must be 1,2 or 3" );
1271 }
1272
1273 // (2) inputFields is (C,F,P,D,D) or (F,P,D,D) and 1 <= (dimension(rank-1), (rank-2)) <= 3 is required.
1274 INTREPID2_TEST_FOR_EXCEPTION( inputFields.rank() < 4 ||
1275 inputFields.rank() > 5, std::invalid_argument,
1276 ">>> ERROR (ArrayTools::matmatProductDataField): inputFields must have rank 4 or 5");
1277 INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(inputFields.rank()-1) < 1 ||
1278 inputFields.extent(inputFields.rank()-1) > 3, std::invalid_argument,
1279 ">>> ERROR (ArrayTools::matmatProductDataField): inputFields dimension (rank-1) must be 1,2 or 3" );
1280 INTREPID2_TEST_FOR_EXCEPTION( inputFields.extent(inputFields.rank()-2) < 1 ||
1281 inputFields.extent(inputFields.rank()-2) > 3, std::invalid_argument,
1282 ">>> ERROR (ArrayTools::matmatProductDataField): inputFields dimension (rank-2) must be 1,2 or 3" );
1283
1284 // (3) outputFields is (C,F,P,D,D)
1285 INTREPID2_TEST_FOR_EXCEPTION( outputFields.rank() != 5, std::invalid_argument,
1286 ">>> ERROR (ArrayTools::matmatProductDataField): outputFields must have rank 5" );
1287 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(3) < 1 ||
1288 outputFields.extent(3) > 3, std::invalid_argument,
1289 ">>> ERROR (ArrayTools::matmatProductDataField): outputFields dimension (3) must be 1,2 or 3" );
1290 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(4) < 1 ||
1291 outputFields.extent(4) > 3, std::invalid_argument,
1292 ">>> ERROR (ArrayTools::matmatProductDataField): outputFields dimension (4) must be 1,2 or 3" );
1293
1294 /*
1295 * Dimension cross-checks:
1296 * (1) inputData vs. inputFields
1297 * (2) outputFields vs. inputData
1298 * (3) outputFields vs. inputFields
1299 *
1300 * Cross-check (2): outputFields(C,F,P,D,D) vs. inputData(C,P), (C,P,D) or (C,P,D,D):
1301 * dimensions C, and D must match in all cases, dimension P must match only when non-constant
1302 * data is specified (P>1). Do not check P dimensions with constant data, i.e., when P=1 in
1303 * inputData(C,1,...)
1304 */
1305 if (inputData.extent(1) > 1) { // check P dimension if P>1 in inputData
1306 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(2) != inputData.extent(1), std::invalid_argument,
1307 ">>> ERROR (ArrayTools::matmatProductDataField): outputFields dimension (2) does not match to inputData dimension (1)" );
1308 }
1309 if (inputData.rank() == 2) { // inputData(C,P) -> C match
1310 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(0) != inputData.extent(0), std::invalid_argument,
1311 ">>> ERROR (ArrayTools::matmatProductDataField): outputFields dimension (0) does not match to inputData dimension (0)" );
1312 }
1313 if (inputData.rank() == 3) { // inputData(C,P,D) -> C, D match
1314 const size_type f1[] = { 0, 3, 4 }, f2[] = { 0, 2, 2 };
1315 for (size_type i=0; i<3; ++i) {
1316 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputData.extent(f2[i]), std::invalid_argument,
1317 ">>> ERROR (ArrayTools::matmatProductDataField): outputFields dimension does not match to inputData dimension" );
1318 }
1319 }
1320 if (inputData.rank() == 4) { // inputData(C,P,D,D) -> C, D, D match
1321 const size_type f1[] = { 0, 3, 4 }, f2[] = { 0, 2, 3 };
1322 for (size_type i=0; i<3; ++i) {
1323 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputData.extent(f2[i]), std::invalid_argument,
1324 ">>> ERROR (ArrayTools::matmatProductDataField): outputFields dimension does not match to inputData dimension" );
1325 }
1326 }
1327
1328 /*
1329 * Cross-checks (1,3):
1330 */
1331 if (inputFields.rank() == 5) {
1332 // 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
1333 if (inputData.extent(1) > 1) { // check P dimension if P>1 in inputData
1334 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(1) != inputFields.extent(2), std::invalid_argument,
1335 ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension (1) does not match to inputFields dimension (2)" );
1336 }
1337 if (inputData.rank() == 2) { // inputData(C,P) -> C match
1338 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(0) != inputFields.extent(0), std::invalid_argument,
1339 ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension (0) does not match to inputFields dimension (0)" );
1340 }
1341 if (inputData.rank() == 3) { // inputData(C,P,D) -> C, D match
1342
1343 const size_type f1[] = { 0, 2, 2 }, f2[] = { 0, 3, 4 };
1344 for (size_type i=0; i<3; ++i) {
1345 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
1346 ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension does not match to inputFields dimension" );
1347 }
1348 }
1349 if (inputData.rank() == 4) { // inputData(C,P,D,D) -> C, D, D match
1350 const size_type f1[] = { 0, 2, 3 }, f2[] = { 0, 3, 4 };
1351 for (size_type i=0; i<3; ++i) {
1352 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
1353 ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension does not match to inputFields dimension" );
1354 }
1355 }
1356
1357 // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(C,F,P,D,D): all dimensions C, F, P, D must match
1358 for (size_type i=0; i<outputFields.rank(); ++i) {
1359 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(i) != inputFields.extent(i), std::invalid_argument,
1360 ">>> ERROR (ArrayTools::matmatProductDataField): outputFields dimension does not match to inputFields dimension" );
1361 }
1362 } else {
1363 // 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
1364 if (inputData.extent(1) > 1) { // check P if P>1 in inputData
1365 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(1) != inputFields.extent(1), std::invalid_argument,
1366 ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension (1) does not match to inputFields dimension (1)" );
1367 }
1368 if (inputData.rank() == 3) {
1369 const size_type f1[] = { 2, 2 }, f2[] = { 2, 3 };
1370 for (size_type i=0; i<2; ++i) {
1371 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
1372 ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension does not match to inputFields dimension" );
1373 }
1374 }
1375 if (inputData.rank() == 4) {
1376 const size_type f1[] = { 2, 3 }, f2[] = { 2, 3 };
1377 for (size_type i=0; i<2; ++i) {
1378 INTREPID2_TEST_FOR_EXCEPTION( inputData.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
1379 ">>> ERROR (ArrayTools::matmatProductDataField): inputData dimension does not match to inputFields dimension" );
1380 }
1381 }
1382
1383 // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(F,P,D,D): dimensions F, P, D must match
1384 {
1385 const size_type f1[] = { 1, 2, 3, 4 }, f2[] = { 0, 1, 2, 3 };
1386 for (size_type i=0; i<4; ++i) {
1387 INTREPID2_TEST_FOR_EXCEPTION( outputFields.extent(f1[i]) != inputFields.extent(f2[i]), std::invalid_argument,
1388 ">>> ERROR (ArrayTools::matmatProductDataField): outputFields dimension does not match to inputFields dimension" );
1389 }
1390 }
1391 }
1392 }
1393#endif
1394 // body
1395 ArrayTools<DeviceType>::Internal::matmatProduct( outputFields,
1396 inputData,
1397 inputFields,
1398 true,
1399 transpose == 't' || transpose == 'T' );
1400 }
1401
1402
1403
1404 template<typename DeviceType>
1405 template<typename outputDataValueType, class ...outputDataProperties,
1406 typename inputDataLeftValueType, class ...inputDataLeftProperties,
1407 typename inputDataRightValueType, class ...inputDataRightProperties>
1408 void
1409 ArrayTools<DeviceType>::matmatProductDataData( Kokkos::DynRankView<outputDataValueType, outputDataProperties...> outputData,
1410 const Kokkos::DynRankView<inputDataLeftValueType, inputDataLeftProperties...> inputDataLeft,
1411 const Kokkos::DynRankView<inputDataRightValueType,inputDataRightProperties...> inputDataRight,
1412 const char transpose ) {
1413
1414#ifdef HAVE_INTREPID2_DEBUG
1415 {
1416 /*
1417 * Check array rank and spatial dimension range (if applicable)
1418 * (1) inputDataLeft(C,P), (C,P,D) or (C,P,D,D); P=1 is admissible to allow multiply by const. left data
1419 * (2) inputDataRight(C,P,D,D) or (P,D,D);
1420 * (3) outputData(C,P,D,D)
1421 */
1422 // (1) inputDataLeft is (C,P), (C,P,D) or (C,P,D,D) and 1 <= D <= 3 is required
1423 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.rank() < 2 ||
1424 inputDataLeft.rank() > 4, std::invalid_argument,
1425 ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft must have rank 2,3 or 4" );
1426 if (inputDataLeft.rank() > 2) {
1427 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(2) < 1 ||
1428 inputDataLeft.extent(2) > 3, std::invalid_argument,
1429 ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft dimension (2) must be 1,2 or 3" );
1430 }
1431 if (inputDataLeft.rank() == 4) {
1432 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(3) < 1 ||
1433 inputDataLeft.extent(3) > 3, std::invalid_argument,
1434 ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft dimension (3) must be 1,2 or 3" );
1435 }
1436
1437 // (2) inputDataRight is (C,P,D,D) or (P,D,D) and 1 <= (D=dimension(rank-1),(rank-2)) <= 3 is required.
1438 INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.rank() < 3 ||
1439 inputDataRight.rank() > 4, std::invalid_argument,
1440 ">>> ERROR (ArrayTools::matmatProductDataData): inputDataRight must have rank 3 or 4" );
1441 INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.extent(inputDataRight.rank()-1) < 1 ||
1442 inputDataRight.extent(inputDataRight.rank()-1) > 3, std::invalid_argument,
1443 ">>> ERROR (ArrayTools::matmatProductDataData): inputDataRight (rank-1) must be 1,2 or 3" );
1444 INTREPID2_TEST_FOR_EXCEPTION( inputDataRight.extent(inputDataRight.rank()-2) < 1 ||
1445 inputDataRight.extent(inputDataRight.rank()-2) > 3, std::invalid_argument,
1446 ">>> ERROR (ArrayTools::matmatProductDataData): inputDataRight (rank-2) must be 1,2 or 3" );
1447
1448 // (3) outputData is (C,P,D,D)
1449 INTREPID2_TEST_FOR_EXCEPTION( outputData.rank() != 4, std::invalid_argument,
1450 ">>> ERROR (ArrayTools::matmatProductDataData): outputData must have rank 4" );
1451 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(2) < 1 ||
1452 outputData.extent(2) > 3, std::invalid_argument,
1453 ">>> ERROR (ArrayTools::matmatProductDataData): outputData (2) must be 1,2 or 3" );
1454 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(3) < 1 ||
1455 outputData.extent(3) > 3, std::invalid_argument,
1456 ">>> ERROR (ArrayTools::matmatProductDataData): outputData (3) must be 1,2 or 3" );
1457
1458 /*
1459 * Dimension cross-checks:
1460 * (1) inputDataLeft vs. inputDataRight
1461 * (2) outputData vs. inputDataLeft
1462 * (3) outputData vs. inputDataRight
1463 *
1464 * Cross-check (2): outputData(C,P,D,D) vs. inputDataLeft(C,P), (C,P,D) or (C,P,D,D):
1465 * dimensions C, and D must match in all cases, dimension P must match only when non-constant
1466 * data is specified (P>1). Do not check P dimensions with constant left data, i.e., when P=1 in
1467 * inputDataLeft(C,1,...)
1468 */
1469 if (inputDataLeft.extent(1) > 1) { // check P dimension if P>1 in inputDataLeft
1470 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(1) != inputDataLeft.extent(1), std::invalid_argument,
1471 ">>> ERROR (ArrayTools::matmatProductDataData): outputData dimension (1) does not match to inputDataLeft dimension (1)" );
1472 }
1473 if (inputDataLeft.rank() == 2) { // inputDataLeft(C,P): check C
1474 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(0) != inputDataLeft.extent(0), std::invalid_argument,
1475 ">>> ERROR (ArrayTools::matmatProductDataData): outputData dimension (0) does not match to inputDataLeft dimension (0)" );
1476 }
1477 if (inputDataLeft.rank() == 3) { // inputDataLeft(C,P,D): check C and D
1478 const size_type f1[] = { 0, 2, 3 }, f2[] = { 0, 2, 2 };
1479 for (size_type i=0; i<3; ++i) {
1480 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataLeft.extent(f2[i]), std::invalid_argument,
1481 ">>> ERROR (ArrayTools::matmatProductDataData): outputData dimension does not match to inputDataLeft dimension" );
1482 }
1483 }
1484 if (inputDataLeft.rank() == 4) { // inputDataLeft(C,P,D,D): check C and D
1485 const size_type f1[] = { 0, 2, 3 }, f2[] = { 0, 2, 3 };
1486 for (size_type i=0; i<3; ++i) {
1487 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataLeft.extent(f2[i]), std::invalid_argument,
1488 ">>> ERROR (ArrayTools::matmatProductDataData): outputData dimension does not match to inputDataLeft dimension" );
1489 }
1490 }
1491
1492 /*
1493 * Cross-checks (1,3):
1494 */
1495 if (inputDataRight.rank() == 4) {
1496 // 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
1497 if (inputDataLeft.extent(1) > 1) { // check P dimension if P>1 in inputDataLeft
1498 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(1) != inputDataRight.extent(1), std::invalid_argument,
1499 ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft dimension (1) does not match to inputDataRight dimension (1)" );
1500 }
1501 if (inputDataLeft.rank() == 2) { // inputDataLeft(C,P): check C
1502 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(0) != inputDataRight.extent(0), std::invalid_argument,
1503 ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft dimension (0) does not match to inputDataRight dimension (0)" );
1504 }
1505 if (inputDataLeft.rank() == 3) { // inputDataLeft(C,P,D): check C and D
1506 const size_type f1[] = { 0, 2, 2 }, f2[] = { 0, 2, 3 };
1507 for (size_type i=0; i<3; ++i) {
1508 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1509 ">>> ERROR (ArrayTools::matmatProductDataData): intputDataLeft dimension does not match to inputDataRight dimension" );
1510 }
1511 }
1512 if (inputDataLeft.rank() == 4) { // inputDataLeft(C,P,D,D): check C and D
1513 const size_type f1[] = { 0, 2, 3 }, f2[] = { 0, 2, 3 };
1514 for (size_type i=0; i<3; ++i) {
1515 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1516 ">>> ERROR (ArrayTools::matmatProductDataData): intputDataLeft dimension does not match to inputDataRight dimension" );
1517 }
1518 }
1519
1520 // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(C,P,D,D): all dimensions C, P, D must match
1521 for (size_type i=0; i< outputData.rank(); ++i) {
1522 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(i) != inputDataRight.extent(i), std::invalid_argument,
1523 ">>> ERROR (ArrayTools::matmatProductDataData): outputData dimension does not match to inputDataRight dimension" );
1524 }
1525 } else {
1526 // Cross-check (1): inputDataLeft(C,P), (C,P,D), or (C,P,D,D) vs. inputDataRight(P,D,D): dimensions P, D must match
1527 if (inputDataLeft.extent(1) > 1) { // check P if P>1 in inputData
1528 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(1) != inputDataRight.extent(0), std::invalid_argument,
1529 ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft dimension (1) does not match to inputDataRight dimension (0)" );
1530 }
1531 if (inputDataLeft.rank() == 3) { // inputDataLeft(C,P,D): check D
1532 const size_type f1[] = { 2, 2 }, f2[] = { 1, 2 };
1533 for (size_type i=0; i<2; ++i) {
1534 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1535 ">>> ERROR (ArrayTools::matmatProductDataData): intputDataLeft dimension does not match to inputDataRight dimension" );
1536 }
1537 }
1538 if (inputDataLeft.rank() == 4) { // inputDataLeft(C,P,D,D): check D
1539 const size_type f1[] = { 2, 3 }, f2[] = { 1, 2 };
1540 for (size_type i=0; i<2; ++i) {
1541 INTREPID2_TEST_FOR_EXCEPTION( inputDataLeft.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1542 ">>> ERROR (ArrayTools::matmatProductDataData): intputDataLeft dimension does not match to inputDataRight dimension" );
1543 }
1544 }
1545 // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(P,D,D): dimensions P, D must match
1546 {
1547 const size_type f1[] = { 1, 2, 3 }, f2[] = { 0, 1, 2 };
1548 for (size_type i=0; i<3; ++i) {
1549 INTREPID2_TEST_FOR_EXCEPTION( outputData.extent(f1[i]) != inputDataRight.extent(f2[i]), std::invalid_argument,
1550 ">>> ERROR (ArrayTools::matmatProductDataData): outputData dimension does not match to inputDataRight dimension" );
1551 }
1552 }
1553 }
1554 }
1555#endif
1556 // body
1557 ArrayTools<DeviceType>::Internal::matmatProduct( outputData,
1558 inputDataLeft,
1559 inputDataRight,
1560 false,
1561 transpose == 't' || transpose == 'T' );
1562 }
1563
1564} // end namespace Intrepid2
1565#endif
Intrepid2::ArrayTools
Utility class that provides methods for higher-order algebraic manipulation of user-defined arrays,...
Definition: Intrepid2_ArrayTools.hpp:77
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:145
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:1245
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:572
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:1409
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:969
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:465
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:815
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:268
Intrepid2::FunctorArrayTools::F_crossProduct
Functor for crossProduct see Intrepid2::ArrayTools for more.
Definition: Intrepid2_ArrayToolsDefTensor.hpp:59
Intrepid2::FunctorArrayTools::F_matmatProduct
Functor for matmatProduct see Intrepid2::ArrayTools for more.
Definition: Intrepid2_ArrayToolsDefTensor.hpp:1123
Intrepid2::FunctorArrayTools::F_matvecProduct
Functor for matvecProduct see Intrepid2::ArrayTools for more.
Definition: Intrepid2_ArrayToolsDefTensor.hpp:677
Intrepid2::FunctorArrayTools::F_outerProduct
Functor for outerProduct see Intrepid2::ArrayTools for more.
Definition: Intrepid2_ArrayToolsDefTensor.hpp:385
Generated by doxygen 1.9.4