Intrepid2
Intrepid2_RealSpaceToolsDef.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_REALSPACETOOLS_DEF_HPP__
17#define __INTREPID2_REALSPACETOOLS_DEF_HPP__
18
19namespace Intrepid2 {
20
21 // ------------------------------------------------------------------------------------
22
23 //
24 // serial version code
25 //
26
27 // ------------------------------------------------------------------------------------
28
29 template<typename DeviceType>
30 template<typename inVecValueType, class ...inVecProperties>
31 KOKKOS_INLINE_FUNCTION
32 inVecValueType
33 RealSpaceTools<DeviceType>::Serial::
34 vectorNorm( const Kokkos::DynRankView<inVecValueType,inVecProperties...> inVec,
35 const ENorm normType ) {
36#ifdef HAVE_INTREPID2_DEBUG
37 {
38 bool dbgInfo = false;
39 INTREPID2_TEST_FOR_DEBUG_ABORT( !(inVec.rank() >= 1 && inVec.rank() <= 5), dbgInfo,
40 ">>> ERROR (RealSpaceTools::vectorNorm): Vector argument must have rank 1!" );
41#ifdef INTREPID2_TEST_FOR_DEBUG_ABORT_OVERRIDE_TO_CONTINUE
42 if (dbgInfo) return inVecValueType(0);
43#endif
44 }
45#endif
46 typedef inVecValueType value_type;
47
48 value_type norm(0);
49
50 const ordinal_type iend = inVec.extent(0);
51 const ordinal_type jend = inVec.extent(1);
52 const ordinal_type kend = inVec.extent(2);
53 const ordinal_type lend = inVec.extent(3);
54 const ordinal_type mend = inVec.extent(4);
55 switch(normType) {
56 case NORM_TWO:{
57 for (ordinal_type i=0;i<iend;++i)
58 for (ordinal_type j=0;j<jend;++j)
59 for (ordinal_type k=0;k<kend;++k)
60 for (ordinal_type l=0;l<lend;++l)
61 for (ordinal_type m=0;m<mend;++m)
62 norm += inVec.access(i,j,k,l,m)*inVec.access(i,j,k,l,m);
63 norm = sqrt(norm);
64 break;
65 }
66 case NORM_INF:{
67 for (ordinal_type i=0;i<iend;++i)
68 for (ordinal_type j=0;j<jend;++j)
69 for (ordinal_type k=0;k<kend;++k)
70 for (ordinal_type l=0;l<lend;++l)
71 for (ordinal_type m=0;m<mend;++m) {
72 const value_type current = Util<value_type>::abs(inVec.access(i,j,k,l,m));
73 norm = (norm < current ? current : norm);
74 }
75 break;
76 }
77 case NORM_ONE:{
78 for (ordinal_type i=0;i<iend;++i)
79 for (ordinal_type j=0;j<jend;++j)
80 for (ordinal_type k=0;k<kend;++k)
81 for (ordinal_type l=0;l<lend;++l)
82 for (ordinal_type m=0;m<mend;++m)
83 norm += Util<value_type>::abs(inVec.access(i,j,k,l,m));
84 break;
85 }
86 default: {
87 INTREPID2_TEST_FOR_ABORT( ( (normType != NORM_TWO) &&
88 (normType != NORM_INF) &&
89 (normType != NORM_ONE) ),
90 ">>> ERROR (RealSpaceTools::vectorNorm): Invalid argument normType.");
91 }
92 }
93 return norm;
94 }
95
96 // ------------------------------------------------------------------------------------
97
98
99 template<typename DeviceType>
100 template<class MatrixViewType>
101 KOKKOS_INLINE_FUNCTION
102 typename MatrixViewType::value_type
103 RealSpaceTools<DeviceType>::Serial::
104 det( const MatrixViewType inMat ) {
105
106 typedef typename decltype(inMat)::non_const_value_type value_type;
107#ifdef HAVE_INTREPID2_DEBUG
108 {
109 bool dbgInfo = false;
110 INTREPID2_TEST_FOR_DEBUG_ABORT( getFunctorRank( inMat ) != 2, dbgInfo,
111 ">>> ERROR (RealSpaceTools::det): Rank of matrix argument must be 2!");
112 INTREPID2_TEST_FOR_DEBUG_ABORT( inMat.extent(0) != inMat.extent(1), dbgInfo,
113 ">>> ERROR (RealSpaceTools::det): Matrix is not square!");
114 INTREPID2_TEST_FOR_DEBUG_ABORT( inMat.extent(0) < 1 || inMat.extent(0) > 3, dbgInfo,
115 ">>> ERROR (RealSpaceTools::det): Spatial dimension must be 1, 2, or 3!");
116#ifdef INTREPID2_TEST_FOR_DEBUG_ABORT_OVERRIDE_TO_CONTINUE
117 if (dbgInfo) return value_type(0);
118#endif
119 }
120#endif
121 const auto dim = inMat.extent(0);
122
123 value_type r_val = 0.0;
124 switch (dim) {
125 case 3:
126 r_val = ( inMat(0,0) * inMat(1,1) * inMat(2,2) +
127 inMat(1,0) * inMat(2,1) * inMat(0,2) +
128 inMat(2,0) * inMat(0,1) * inMat(1,2) -
129 inMat(2,0) * inMat(1,1) * inMat(0,2) -
130 inMat(0,0) * inMat(2,1) * inMat(1,2) -
131 inMat(1,0) * inMat(0,1) * inMat(2,2) );
132 break;
133 case 2:
134 r_val = ( inMat(0,0) * inMat(1,1) -
135 inMat(0,1) * inMat(1,0) );
136 break;
137 case 1:
138 r_val = ( inMat(0,0) );
139 break;
140 }
141 return r_val;
142 }
143
144 // ------------------------------------------------------------------------------------
145
146 template<typename DeviceType>
147 template<typename inVec1ValueType, class ...inVec1Properties,
148 typename inVec2ValueType, class ...inVec2Properties>
149 KOKKOS_INLINE_FUNCTION
150 inVec1ValueType
151 RealSpaceTools<DeviceType>::Serial::
152 dot( const Kokkos::DynRankView<inVec1ValueType,inVec1Properties...> inVec1,
153 const Kokkos::DynRankView<inVec2ValueType,inVec2Properties...> inVec2 ) {
154
155#ifdef HAVE_INTREPID2_DEBUG
156 {
157 bool dbgInfo = false;
158 INTREPID2_TEST_FOR_DEBUG_ABORT( inVec1.rank() != 1 || inVec2.rank() != 1, dbgInfo,
159 ">>> ERROR (RealSpaceTools::dot): Vector arguments must have rank 1!");
160 INTREPID2_TEST_FOR_DEBUG_ABORT( inVec1.extent(0) != inVec2.extent(0), dbgInfo,
161 ">>> ERROR (RealSpaceTools::dot): Dimensions of vector arguments must agree!");
162#ifdef INTREPID2_TEST_FOR_EXCEPTION_OVERRIDE_TO_CONTINUE
163 if (dbgInfo) return inVec1ValueType(0);
164#endif
165 }
166#endif
167 typedef inVec1ValueType value_type;
168
169 // designed for small problems
170 value_type r_val(0);
171
172 const ordinal_type iend = inVec1.extent(0);
173 for (ordinal_type i=0;i<iend;++i)
174 r_val += inVec1(i)*inVec2(i);
175
176 return r_val;
177 }
178
179 // ------------------------------------------------------------------------------------
180
181 //
182 // use parallel for
183 //
184
185 // ------------------------------------------------------------------------------------
186
187 namespace FunctorRealSpaceTools {
188
192 template<typename OutputViewType,
193 typename inputViewType>
194 struct F_extractScalarValues {
195 OutputViewType _output;
196 inputViewType _input;
197
198 KOKKOS_INLINE_FUNCTION
199 F_extractScalarValues( OutputViewType output_,
200 inputViewType input_ )
201 : _output(output_), _input(input_) {}
202
203 KOKKOS_INLINE_FUNCTION
204 void operator()(const ordinal_type i) const {
205 const ordinal_type jend = _input.extent(1);
206 const ordinal_type kend = _input.extent(2);
207 const ordinal_type lend = _input.extent(3);
208 const ordinal_type mend = _input.extent(4);
209
210 for (ordinal_type j=0;j<jend;++j)
211 for (ordinal_type k=0;k<kend;++k)
212 for (ordinal_type l=0;l<lend;++l)
213 for (ordinal_type m=0;m<mend;++m)
214 _output.access(i,j,k,l,m) = get_scalar_value(_input.access(i,j,k,l,m));
215 }
216 };
217 }
218
219 template<typename DeviceType>
220 template<typename outputValueType, class ...outputProperties,
221 typename inputValueType, class ...inputProperties>
222 void
223 RealSpaceTools<DeviceType>::
224 extractScalarValues( Kokkos::DynRankView<outputValueType,outputProperties...> output,
225 const Kokkos::DynRankView<inputValueType, inputProperties...> input ) {
226
227 using MemSpaceType = typename DeviceType::memory_space;
228 constexpr bool are_accessible =
229 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
230 typename decltype(output)::memory_space>::accessible &&
231 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
232 typename decltype(input)::memory_space>::accessible;
233 static_assert(are_accessible, "RealSpaceTools<DeviceType>::extractScalarValues(..): input/output views' memory spaces are not compatible with DeviceType");
234
235 using FunctorType = FunctorRealSpaceTools::F_extractScalarValues<decltype(output),decltype(input)>;
236
237 const auto loopSize = input.extent(0);
238 using ExecSpaceType = typename DeviceType::execution_space;
239 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
240 Kokkos::parallel_for( policy, FunctorType(output, input) );
241 }
242
243 namespace FunctorRealSpaceTools {
247 template<typename OutputViewType,
248 typename inputViewType>
249 struct F_clone {
250 OutputViewType _output;
251 inputViewType _input;
252
253 KOKKOS_INLINE_FUNCTION
254 F_clone( OutputViewType output_,
255 inputViewType input_ )
256 : _output(output_), _input(input_) {}
257
258 KOKKOS_INLINE_FUNCTION
259 void operator()(const ordinal_type iter) const {
260 ordinal_type k0(0), k1(0), k2(0);
261 const auto rankDiff = _output.rank() - _input.rank();
262 switch (rankDiff) {
263 case 3:
264 unrollIndex( k0, k1, k2,
265 _output.extent(0),
266 _output.extent(1),
267 _output.extent(2),
268 iter );
269 break;
270 case 2:
271 unrollIndex( k0, k1,
272 _output.extent(0),
273 _output.extent(1),
274 iter );
275 break;
276 case 1:
277 k0 = iter;
278 break;
279 }
280
281 auto out = (rankDiff == 3 ? Kokkos::subview(_output, k0, k1, k2, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()) :
282 rankDiff == 2 ? Kokkos::subview(_output, k0, k1, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()) :
283 rankDiff == 1 ? Kokkos::subview(_output, k0, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()) :
284 Kokkos::subview(_output, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()));
285 const ordinal_type iend = _input.extent(0);
286 const ordinal_type jend = _input.extent(1);
287 const ordinal_type kend = _input.extent(2);
288 for (ordinal_type i=0;i<iend;++i)
289 for (ordinal_type j=0;j<jend;++j)
290 for (ordinal_type k=0;k<kend;++k)
291 out.access(i,j,k) = _input.access(i,j,k);
292 }
293 };
294 }
295
296 template<typename DeviceType>
297 template<typename outputValueType, class ...outputProperties,
298 typename inputValueType, class ...inputProperties>
299 void
300 RealSpaceTools<DeviceType>::
301 clone( Kokkos::DynRankView<outputValueType,outputProperties...> output,
302 const Kokkos::DynRankView<inputValueType,inputProperties...> input ) {
303#ifdef HAVE_INTREPID2_DEBUG
304 {
305 // input has at most rank 3
306 INTREPID2_TEST_FOR_EXCEPTION( input.rank() > 3, std::invalid_argument,
307 ">>> ERROR (RealSpaceTools::clone): Input container has rank larger than 3.");
308
309
310 INTREPID2_TEST_FOR_EXCEPTION( input.rank() > output.rank(), std::invalid_argument,
311 ">>> ERROR (RealSpaceTools::clone): Input container rank should be smaller than ouput rank.");
312
313 const size_type rankDiff = output.rank() - input.rank();
314
315 // Difference of output and input rank is at most 3.
316 INTREPID2_TEST_FOR_EXCEPTION( rankDiff > 3, std::invalid_argument,
317 ">>> ERROR (RealSpaceTools::clone): Difference between output container rank and input container rank is larger than 3.");
318
319
320 for (size_type i=0;i<input.rank();++i) {
321 INTREPID2_TEST_FOR_EXCEPTION( input.extent(i) != output.extent(rankDiff + i), std::invalid_argument,
322 ">>> ERROR (RealSpaceTools::clone): Dimensions of array arguments do not agree!");
323 }
324 }
325#endif
326 using MemSpaceType = typename DeviceType::memory_space;
327 constexpr bool are_accessible =
328 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
329 typename decltype(output)::memory_space>::accessible &&
330 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
331 typename decltype(input)::memory_space>::accessible;
332 static_assert(are_accessible, "RealSpaceTools<DeviceType>::clone(..): input/output views' memory spaces are not compatible with DeviceType");
333
334 using FunctorType = FunctorRealSpaceTools::F_clone<decltype(output),decltype(input)>;
335
336 size_type loopSize = 1;
337 const ordinal_type out_rank = output.rank();
338 const ordinal_type in_rank = input.rank();
339 const ordinal_type rankDiff = out_rank - in_rank;
340 for (ordinal_type i=0;i<rankDiff;++i)
341 loopSize *= output.extent(i);
342
343 using ExecSpaceType = typename DeviceType::execution_space;
344 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
345 Kokkos::parallel_for( policy, FunctorType(output, input) );
346 }
347
348 namespace FunctorRealSpaceTools {
352 template<typename absArrayViewType,
353 typename inArrayViewType>
354 struct F_absval {
355 absArrayViewType _absArray;
356 inArrayViewType _inArray;
357
358 KOKKOS_INLINE_FUNCTION
359 F_absval( absArrayViewType absArray_,
360 inArrayViewType inArray_ )
361 : _absArray(absArray_), _inArray(inArray_) {}
362
363 KOKKOS_INLINE_FUNCTION
364 void operator()(const ordinal_type i) const {
365 const ordinal_type jend = _inArray.extent(1);
366 const ordinal_type kend = _inArray.extent(2);
367 const ordinal_type lend = _inArray.extent(3);
368 const ordinal_type mend = _inArray.extent(4);
369
370 typedef typename inArrayViewType::non_const_value_type value_type;
371
372 for (ordinal_type j=0;j<jend;++j)
373 for (ordinal_type k=0;k<kend;++k)
374 for (ordinal_type l=0;l<lend;++l)
375 for (ordinal_type m=0;m<mend;++m)
376 _absArray.access(i,j,k,l,m) = Util<value_type>::abs(_inArray.access(i,j,k,l,m));
377 }
378 };
379 }
380
381 template<typename DeviceType>
382 template<typename absArrayValueType, class ...absArrayProperties,
383 typename inArrayValueType, class ...inArrayProperties>
384 void
385 RealSpaceTools<DeviceType>::
386 absval( Kokkos::DynRankView<absArrayValueType,absArrayProperties...> absArray,
387 const Kokkos::DynRankView<inArrayValueType, inArrayProperties...> inArray ) {
388#ifdef HAVE_INTREPID2_DEBUG
389 {
390 INTREPID2_TEST_FOR_EXCEPTION( inArray.rank() > 5, std::invalid_argument,
391 ">>> ERROR (RealSpaceTools::absval): Input array container has rank larger than 5.");
392
393 INTREPID2_TEST_FOR_EXCEPTION( inArray.rank() != absArray.rank(), std::invalid_argument,
394 ">>> ERROR (RealSpaceTools::absval): Array arguments must have identical ranks!");
395 for (size_type i=0;i<inArray.rank();++i) {
396 INTREPID2_TEST_FOR_EXCEPTION( inArray.extent(i) != absArray.extent(i), std::invalid_argument,
397 ">>> ERROR (RealSpaceTools::absval): Dimensions of array arguments do not agree!");
398 }
399 }
400#endif
401 using MemSpaceType = typename DeviceType::memory_space;
402 constexpr bool are_accessible =
403 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
404 typename decltype(absArray)::memory_space>::accessible &&
405 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
406 typename decltype(inArray)::memory_space>::accessible;
407 static_assert(are_accessible, "RealSpaceTools<DeviceType>::absval(..): input/output views' memory spaces are not compatible with DeviceType");
408
409 using FunctorType = FunctorRealSpaceTools::F_absval<decltype(absArray),decltype(inArray)>;
410
411 const auto loopSize = inArray.extent(0);
412 using ExecSpaceType = typename DeviceType::execution_space;
413 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
414 Kokkos::parallel_for( policy, FunctorType(absArray, inArray) );
415 }
416
417 // ------------------------------------------------------------------------------------
418
419 template<typename DeviceType>
420 template<typename inoutArrayValueType, class ...inoutAbsArrayProperties>
421 void
422 RealSpaceTools<DeviceType>::
423 absval( Kokkos::DynRankView<inoutArrayValueType,inoutAbsArrayProperties...> inoutAbsArray ) {
424 RealSpaceTools<DeviceType>::absval(inoutAbsArray, inoutAbsArray);
425 }
426
427 // ------------------------------------------------------------------------------------
428
429 namespace FunctorRealSpaceTools {
433 template<typename normArrayViewType,
434 typename inVecViewType>
435 struct F_vectorNorm {
436 normArrayViewType _normArray;
437 inVecViewType _inVecs;
438 const ENorm _normType;
439
440 KOKKOS_INLINE_FUNCTION
441 F_vectorNorm( normArrayViewType normArray_,
442 inVecViewType inVecs_,
443 const ENorm normType_ )
444 : _normArray(normArray_), _inVecs(inVecs_), _normType(normType_) {}
445
446 KOKKOS_INLINE_FUNCTION
447 void operator()(const ordinal_type iter) const {
448 ordinal_type i, j;
449 unrollIndex( i, j,
450 _normArray.extent(0),
451 _normArray.extent(1),
452 iter );
453
454 auto vec = ( _inVecs.rank() == 2 ? Kokkos::subview(_inVecs, i, Kokkos::ALL()) :
455 Kokkos::subview(_inVecs, i, j, Kokkos::ALL()) );
456
457 _normArray.access(i, j) = RealSpaceTools<>::Serial::vectorNorm(vec, _normType);
458 }
459 };
460 }
461
462 template<typename DeviceType>
463 template<typename normArrayValueType, class ...normArrayProperties,
464 typename inVecValueType, class ...inVecProperties>
465 void
466 RealSpaceTools<DeviceType>::
467 vectorNorm( Kokkos::DynRankView<normArrayValueType,normArrayProperties...> normArray,
468 const Kokkos::DynRankView<inVecValueType, inVecProperties...> inVecs,
469 const ENorm normType ) {
470#ifdef HAVE_INTREPID2_DEBUG
471 {
472 INTREPID2_TEST_FOR_EXCEPTION( inVecs.rank() != (normArray.rank()+1), std::invalid_argument,
473 ">>> ERROR (RealSpaceTools::vectorNorm): Ranks of norm and vector array arguments are incompatible!");
474 INTREPID2_TEST_FOR_EXCEPTION( inVecs.rank() < 2 || inVecs.rank() > 3, std::invalid_argument,
475 ">>> ERROR (RealSpaceTools::vectorNorm): Rank of vector array must be 2 or 3!");
476 for (size_type i=0;i<inVecs.rank()-1;++i) {
477 INTREPID2_TEST_FOR_EXCEPTION( inVecs.extent(i) != normArray.extent(i), std::invalid_argument,
478 ">>> ERROR (RealSpaceTools::vectorNorm): Dimensions of norm and vector arguments do not agree!");
479 }
480 }
481#endif
482
483 using MemSpaceType = typename DeviceType::memory_space;
484 constexpr bool are_accessible =
485 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
486 typename decltype(normArray)::memory_space>::accessible &&
487 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
488 typename decltype(inVecs)::memory_space>::accessible;
489 static_assert(are_accessible, "RealSpaceTools<DeviceType>::vectorNorm(..): input/output views' memory spaces are not compatible with DeviceType");
490 using FunctorType = FunctorRealSpaceTools::F_vectorNorm<decltype(normArray),decltype(inVecs)>;
491
492 // normArray rank is either 1 or 2
493 const auto loopSize = normArray.extent(0)*normArray.extent(1);
494 using ExecSpaceType = typename DeviceType::execution_space;
495 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
496 Kokkos::parallel_for( policy, FunctorType(normArray, inVecs, normType) );
497 }
498
499 // ------------------------------------------------------------------------------------
500
501 namespace FunctorRealSpaceTools {
505 template<typename transposeMatViewType,
506 typename inMatViewType>
507 struct F_transpose {
508 transposeMatViewType _transposeMats;
509 inMatViewType _inMats;
510
511 KOKKOS_INLINE_FUNCTION
512 F_transpose( transposeMatViewType transposeMats_,
513 inMatViewType inMats_)
514 : _transposeMats(transposeMats_), _inMats(inMats_) {}
515
516 KOKKOS_INLINE_FUNCTION
517 void operator()(const size_type iter) const {
518 // the rank of normArray is either 1 or 2
519 const auto r = _transposeMats.rank();
520 ordinal_type _i = iter, _j = 0;
521
522 if ( r > 3 )
523 unrollIndex( _i, _j,
524 _transposeMats.extent(0),
525 _transposeMats.extent(1),
526 iter );
527
528 auto dst = ( r == 2 ? Kokkos::subview(_transposeMats, Kokkos::ALL(), Kokkos::ALL()) :
529 r == 3 ? Kokkos::subview(_transposeMats, _i, Kokkos::ALL(), Kokkos::ALL()) :
530 Kokkos::subview(_transposeMats, _i, _j, Kokkos::ALL(), Kokkos::ALL()) );
531
532 auto src = ( r == 2 ? Kokkos::subview(_inMats, Kokkos::ALL(), Kokkos::ALL()) :
533 r == 3 ? Kokkos::subview(_inMats, _i, Kokkos::ALL(), Kokkos::ALL()) :
534 Kokkos::subview(_inMats, _i, _j, Kokkos::ALL(), Kokkos::ALL()) );
535
536 for (size_type i=0;i<src.extent(0);++i) {
537 dst(i, i) = src(i, i);
538 for (size_type j=i+1;j<src.extent(1);++j) {
539 dst(i, j) = src(j, i);
540 dst(j, i) = src(i, j);
541 }
542 }
543 }
544 };
545 }
546
547 template<typename DeviceType>
548 template<typename transposeMatValueType, class ...transposeMatProperties,
549 typename inMatValueType, class ...inMatProperties>
550 void
551 RealSpaceTools<DeviceType>::
552 transpose( Kokkos::DynRankView<transposeMatValueType,transposeMatProperties...> transposeMats,
553 const Kokkos::DynRankView<inMatValueType, inMatProperties...> inMats ) {
554
555#ifdef HAVE_INTREPID2_DEBUG
556 {
557
558 INTREPID2_TEST_FOR_EXCEPTION( inMats.rank() != transposeMats.rank(), std::invalid_argument,
559 ">>> ERROR (RealSpaceTools::transpose): Matrix array arguments do not have identical ranks!");
560 INTREPID2_TEST_FOR_EXCEPTION( inMats.rank() < 2 || inMats.rank() > 4, std::invalid_argument,
561 ">>> ERROR (RealSpaceTools::transpose): Rank of matrix array must be 2, 3, or 4!");
562 for (size_type i=0;i<inMats.rank();++i) {
563 INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(i) != transposeMats.extent(i), std::invalid_argument,
564 ">>> ERROR (RealSpaceTools::transpose): Dimensions of matrix arguments do not agree!");
565 }
566 INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(inMats.rank()-2) != inMats.extent(inMats.rank()-1), std::invalid_argument,
567 ">>> ERROR (RealSpaceTools::transpose): Matrices are not square!");
568 }
569#endif
570
571 using MemSpaceType = typename DeviceType::memory_space;
572 constexpr bool are_accessible =
573 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
574 typename decltype(transposeMats)::memory_space>::accessible &&
575 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
576 typename decltype(inMats)::memory_space>::accessible;
577 static_assert(are_accessible, "RealSpaceTools<DeviceType>::transpose(..): input/output views' memory spaces are not compatible with DeviceType");
578 using FunctorType = FunctorRealSpaceTools::F_transpose<decltype(transposeMats),decltype(inMats)>;
579
580 const auto r = transposeMats.rank();
581 const auto loopSize = ( r == 2 ? 1 :
582 r == 3 ? transposeMats.extent(0) :
583 transposeMats.extent(0)*transposeMats.extent(1) );
584
585 using ExecSpaceType = typename DeviceType::execution_space;
586 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
587 Kokkos::parallel_for( policy, FunctorType(transposeMats, inMats) );
588 }
589
590 // ------------------------------------------------------------------------------------
591
592 namespace FunctorRealSpaceTools {
596 template<typename inverseMatViewType,
597 typename inMatViewType>
598 struct F_inverse {
599 typedef typename inMatViewType::non_const_value_type value_type;
600 inverseMatViewType _inverseMats;
601 inMatViewType _inMats;
602
603 KOKKOS_INLINE_FUNCTION
604 F_inverse( inverseMatViewType inverseMats_,
605 inMatViewType inMats_ )
606 : _inverseMats(inverseMats_), _inMats(inMats_) {}
607
608 template<typename matViewType,
609 typename invViewType>
610 KOKKOS_FORCEINLINE_FUNCTION
611 static void
612 apply_inverse( invViewType inv,
613 const matViewType mat ) {
614 // compute determinant
615 const value_type val = RealSpaceTools<>::Serial::det(mat);
616
617#ifdef HAVE_INTREPID2_DEBUG
618 {
619#ifdef HAVE_INTREPID2_DEBUG_INF_CHECK
620 bool dbgInfo = false;
621 INTREPID2_TEST_FOR_DEBUG_ABORT( val == 0, dbgInfo,
622 ">>> ERROR (Matrix): Inverse of a singular matrix is undefined!");
623#ifdef INTREPID2_TEST_FOR_DEBUG_ABORT_OVERRIDE_TO_CONTINUE
624 if (dbgInfo) return;
625#endif
626#endif
627 }
628#endif
629 switch (mat.extent(0)) {
630 case 1: {
631 inv(0,0) = value_type(1)/mat(0,0);
632 break;
633 }
634 case 2: {
635 inv(0,0) = mat(1,1)/val;
636 inv(1,1) = mat(0,0)/val;
637
638 inv(1,0) = - mat(1,0)/val;
639 inv(0,1) = - mat(0,1)/val;
640 break;
641 }
642 case 3: {
643 value_type val0, val1, val2;
644
645 val0 = mat(1,1)*mat(2,2) - mat(2,1)*mat(1,2);
646 val1 = - mat(1,0)*mat(2,2) + mat(2,0)*mat(1,2);
647 val2 = mat(1,0)*mat(2,1) - mat(2,0)*mat(1,1);
648
649 inv(0,0) = val0/val;
650 inv(1,0) = val1/val;
651 inv(2,0) = val2/val;
652
653 val0 = mat(2,1)*mat(0,2) - mat(0,1)*mat(2,2);
654 val1 = mat(0,0)*mat(2,2) - mat(2,0)*mat(0,2);
655 val2 = - mat(0,0)*mat(2,1) + mat(2,0)*mat(0,1);
656
657 inv(0,1) = val0/val;
658 inv(1,1) = val1/val;
659 inv(2,1) = val2/val;
660
661 val0 = mat(0,1)*mat(1,2) - mat(1,1)*mat(0,2);
662 val1 = - mat(0,0)*mat(1,2) + mat(1,0)*mat(0,2);
663 val2 = mat(0,0)*mat(1,1) - mat(1,0)*mat(0,1);
664
665 inv(0,2) = val0/val;
666 inv(1,2) = val1/val;
667 inv(2,2) = val2/val;
668 break;
669 }
670 }
671 }
672
673 template< bool B, class T = void >
674 using enable_if_t = typename std::enable_if<B,T>::type;
675
676 template<int M=0>
677 KOKKOS_INLINE_FUNCTION
678 enable_if_t<M==0 && supports_rank_4<inMatViewType>::value >
679 operator()(const ordinal_type cl,
680 const ordinal_type pt) const {
681 const auto mat = Kokkos::subview(_inMats, cl, pt, Kokkos::ALL(), Kokkos::ALL());
682 auto inv = Kokkos::subview(_inverseMats, cl, pt, Kokkos::ALL(), Kokkos::ALL());
683
684 apply_inverse( inv, mat );
685 }
686
687 template<int M=0>
688 KOKKOS_INLINE_FUNCTION
689 enable_if_t<M==0 && !supports_rank_4<inMatViewType>::value >
690 operator()(const ordinal_type cl, const ordinal_type pt) const {
691 // empty implementation to allow compilation of parallel_for. (this combination never gets invoked at run-time, but does occur at compile-time.)
692 }
693
694 template<int M=0>
695 KOKKOS_INLINE_FUNCTION
696 enable_if_t<M==0 && supports_rank_3<inMatViewType>::value >
697 operator()(const ordinal_type pt) const {
698 const auto mat = Kokkos::subview(_inMats, pt, Kokkos::ALL(), Kokkos::ALL());
699 auto inv = Kokkos::subview(_inverseMats, pt, Kokkos::ALL(), Kokkos::ALL());
700
701 apply_inverse( inv, mat );
702 }
703
704 template<int M=0>
705 KOKKOS_INLINE_FUNCTION
706 enable_if_t<M==0 && !supports_rank_3<inMatViewType>::value >
707 operator()(const ordinal_type pt) const {
708 // empty implementation to allow compilation of parallel_for. (this combination never gets invoked at run-time, but does occur at compile-time.)
709 }
710 };
711 }
712
713 template<typename DeviceType>
714 template<class InverseMatrixViewType, class MatrixViewType>
715 void
716 RealSpaceTools<DeviceType>::
717 inverse( InverseMatrixViewType inverseMats, MatrixViewType inMats ) {
718 const unsigned rank = getFunctorRank(inMats);
719#ifdef HAVE_INTREPID2_DEBUG
720 {
721 INTREPID2_TEST_FOR_EXCEPTION( rank != getFunctorRank(inverseMats), std::invalid_argument,
722 ">>> ERROR (RealSpaceTools::inverse): Matrix array arguments do not have identical ranks!");
723 INTREPID2_TEST_FOR_EXCEPTION( rank < 3 || rank > 4, std::invalid_argument,
724 ">>> ERROR (RealSpaceTools::inverse): Rank of matrix array must be 3, or 4!");
725 for (size_type i=0;i<rank;++i) {
726 INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(i) != inverseMats.extent(i), std::invalid_argument,
727 ">>> ERROR (RealSpaceTools::inverse): Dimensions of matrix arguments do not agree!");
728 }
729 INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(rank-2) != inMats.extent(rank-1), std::invalid_argument,
730 ">>> ERROR (RealSpaceTools::inverse): Matrices are not square!");
731 INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(rank-2) < 1 || inMats.extent(rank-2) > 3, std::invalid_argument,
732 ">>> ERROR (RealSpaceTools::inverse): Spatial dimension must be 1, 2, or 3!");
733 }
734#endif
735
736 using ExecSpaceType = typename DeviceType::execution_space;
737 using FunctorType = FunctorRealSpaceTools::F_inverse<InverseMatrixViewType, MatrixViewType>;
738
739 switch (rank) {
740 case 3: { // output P,D,D and input P,D,D
741 using range_policy_type = Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> >;
742 range_policy_type policy(0, inverseMats.extent(0));
743 Kokkos::parallel_for( policy, FunctorType(inverseMats, inMats) );
744 break;
745 }
746 case 4: { // output C,P,D,D and input C,P,D,D
747 using range_policy_type = Kokkos::MDRangePolicy
748 < ExecSpaceType, Kokkos::Rank<2>, Kokkos::IndexType<ordinal_type> >;
749 range_policy_type policy( { 0, 0 },
750 { inverseMats.extent(0), inverseMats.extent(1) } );
751 Kokkos::parallel_for( policy, FunctorType(inverseMats, inMats) );
752 break;
753 }
754 default: {
755 INTREPID2_TEST_FOR_EXCEPTION( true, std::invalid_argument,
756 ">>> ERROR (RealSpaceTools::inverse): Rank of matrix array must be 2, 3, or 4!");
757 }
758 }
759 }
760
761 // ------------------------------------------------------------------------------------
762
763 namespace FunctorRealSpaceTools {
767 template<typename detArrayViewType,
768 typename inMatViewType, int rank>
769 struct F_det {
770 detArrayViewType _detArray;
771 inMatViewType _inMats;
772
773 KOKKOS_INLINE_FUNCTION
774 F_det( detArrayViewType detArray_,
775 inMatViewType inMats_ )
776 : _detArray(detArray_), _inMats(inMats_) {}
777
778 template< bool B, class T = void >
779 using enable_if_t = typename std::enable_if<B,T>::type;
780
781 template<int M=rank>
782 KOKKOS_INLINE_FUNCTION
783 enable_if_t<M==1 && supports_rank_3<inMatViewType>::value >
784 operator()(const ordinal_type pt) const {
785 const auto mat = Kokkos::subview(_inMats, pt, Kokkos::ALL(), Kokkos::ALL());
786 _detArray(pt) = RealSpaceTools<>::Serial::det(mat);
787 }
788
789 template<int M=rank>
790 KOKKOS_INLINE_FUNCTION
791 enable_if_t<M==1 && !supports_rank_3<inMatViewType>::value >
792 operator()(const ordinal_type pt) const {
793 // empty implementation to allow compilation of parallel_for. (this combination never gets invoked at run-time, but does occur at compile-time.)
794 }
795
796 template<int M=rank>
797 KOKKOS_INLINE_FUNCTION
798 enable_if_t<M==2 && supports_rank_4<inMatViewType>::value >
799 operator()(const ordinal_type cl,
800 const ordinal_type pt) const {
801 const auto mat = Kokkos::subview(_inMats, cl, pt, Kokkos::ALL(), Kokkos::ALL());
802 _detArray(cl, pt) = RealSpaceTools<>::Serial::det(mat);
803 }
804
805 template<int M=rank>
806 KOKKOS_INLINE_FUNCTION
807 enable_if_t<M==2 && !supports_rank_4<inMatViewType>::value >
808 operator()(const ordinal_type cl,
809 const ordinal_type pt) const {
810 // empty implementation to allow compilation of parallel_for. (this combination never gets invoked at run-time, but does occur at compile-time.)
811 }
812 };
813 }
814
815template<class DeterminantArrayViewType, class MatrixViewType>
816static void
817det( DeterminantArrayViewType detArray, const MatrixViewType inMats );
818
819 template<typename DeviceType>
820 template<class DeterminantArrayViewType, class MatrixViewType>
821 void
822 RealSpaceTools<DeviceType>::
823 det( DeterminantArrayViewType detArray, const MatrixViewType inMats ) {
824
825#ifdef HAVE_INTREPID2_DEBUG
826 {
827 const unsigned matrixRank = getFunctorRank(inMats);
828 const unsigned detRank = getFunctorRank(detArray);
829 INTREPID2_TEST_FOR_EXCEPTION( matrixRank != detRank+2, std::invalid_argument,
830 ">>> ERROR (RealSpaceTools::det): Determinant and matrix array arguments do not have compatible ranks!");
831 INTREPID2_TEST_FOR_EXCEPTION( matrixRank < 3 || matrixRank > 4, std::invalid_argument,
832 ">>> ERROR (RealSpaceTools::det): Rank of matrix array must be 3 or 4!");
833 for (size_type i=0;i<matrixRank-2;++i) {
834 INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(i) != detArray.extent(i), std::invalid_argument,
835 ">>> ERROR (RealSpaceTools::det): Dimensions of determinant and matrix array arguments do not agree!");
836 }
837 INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(matrixRank-2) != inMats.extent(matrixRank-1), std::invalid_argument,
838 ">>> ERROR (RealSpaceTools::det): Matrices are not square!");
839 INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(matrixRank-2) < 1 || inMats.extent(matrixRank-2) > 3, std::invalid_argument,
840 ">>> ERROR (RealSpaceTools::det): Spatial dimension must be 1, 2, or 3!");
841 }
842#endif
843
844 typedef typename DeviceType::execution_space ExecSpaceType;
845
846 const int detArrayRank = getFunctorRank(detArray);
847
848 switch (detArrayRank) {
849 case 1: { // output P and input P,D,D
850 using FunctorType = FunctorRealSpaceTools::F_det<DeterminantArrayViewType,MatrixViewType,1>;
851 using range_policy_type = Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> >;
852 range_policy_type policy(0, detArray.extent(0));
853 Kokkos::parallel_for( policy, FunctorType(detArray, inMats) );
854 break;
855 }
856 case 2: { // output C,P and input C,P,D,D
857 using FunctorType = FunctorRealSpaceTools::F_det<DeterminantArrayViewType,MatrixViewType,2>;
858 using range_policy_type = Kokkos::MDRangePolicy
859 < ExecSpaceType, Kokkos::Rank<2>, Kokkos::IndexType<ordinal_type> >;
860 range_policy_type policy( { 0, 0 },
861 { detArray.extent(0), detArray.extent(1) } );
862 Kokkos::parallel_for( policy, FunctorType(detArray, inMats) );
863 break;
864 }
865 default: {
866 INTREPID2_TEST_FOR_EXCEPTION( true, std::invalid_argument,
867 ">>> ERROR (RealSpaceTools::det): Rank of detArray must be 1 or 2!");
868 }
869 }
870 }
871
872 // ------------------------------------------------------------------------------------
873
874 namespace FunctorRealSpaceTools {
878 template<typename sumArrayViewType,
879 typename inArray1Viewtype,
880 typename inArray2ViewType>
881 struct F_add {
882 sumArrayViewType _sumArray;
883 inArray1Viewtype _inArray1;
884 inArray2ViewType _inArray2;
885
886 KOKKOS_INLINE_FUNCTION
887 F_add( sumArrayViewType sumArray_,
888 inArray1Viewtype inArray1_,
889 inArray2ViewType inArray2_ )
890 : _sumArray(sumArray_), _inArray1(inArray1_), _inArray2(inArray2_) {}
891
892 KOKKOS_INLINE_FUNCTION
893 void operator()(const ordinal_type i) const {
894 const ordinal_type jend = _sumArray.extent(1);
895 const ordinal_type kend = _sumArray.extent(2);
896 const ordinal_type lend = _sumArray.extent(3);
897 const ordinal_type mend = _sumArray.extent(4);
898
899 for (ordinal_type j=0;j<jend;++j)
900 for (ordinal_type k=0;k<kend;++k)
901 for (ordinal_type l=0;l<lend;++l)
902 for (ordinal_type m=0;m<mend;++m)
903 _sumArray.access(i,j,k,l,m) = _inArray1.access(i,j,k,l,m) + _inArray2.access(i,j,k,l,m);
904 }
905 };
906 }
907
908 template<typename DeviceType>
909 template<typename sumArrayValueType, class ...sumArrayProperties,
910 typename inArray1ValueType, class ...inArray1Properties,
911 typename inArray2ValueType, class ...inArray2Properties>
912 void
913 RealSpaceTools<DeviceType>::
914 add( Kokkos::DynRankView<sumArrayValueType,sumArrayProperties...> sumArray,
915 const Kokkos::DynRankView<inArray1ValueType,inArray1Properties...> inArray1,
916 const Kokkos::DynRankView<inArray2ValueType,inArray2Properties...> inArray2 ) {
917
918#ifdef HAVE_INTREPID2_DEBUG
919 {
920 INTREPID2_TEST_FOR_EXCEPTION( inArray1.rank() != inArray2.rank() ||
921 inArray1.rank() != sumArray.rank(), std::invalid_argument,
922 ">>> ERROR (RealSpaceTools::add): Array arguments must have identical ranks!");
923 for (size_type i=0;i<inArray1.rank();++i) {
924 INTREPID2_TEST_FOR_EXCEPTION( inArray1.extent(i) != inArray2.extent(i) ||
925 inArray1.extent(i) != sumArray.extent(i), std::invalid_argument,
926 ">>> ERROR (RealSpaceTools::add): Dimensions of array arguments do not agree!");
927 }
928 }
929#endif
930 using MemSpaceType = typename DeviceType::memory_space;
931 constexpr bool are_accessible =
932 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
933 typename decltype(sumArray)::memory_space>::accessible &&
934 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
935 typename decltype(inArray1)::memory_space>::accessible &&
936 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
937 typename decltype(inArray2)::memory_space>::accessible;
938 static_assert(are_accessible, "RealSpaceTools<DeviceType>::add(..): input/output views' memory spaces are not compatible with DeviceType");
939
940 using FunctorType = FunctorRealSpaceTools::F_add<decltype(sumArray),decltype(inArray1),decltype(inArray2)>;
941
942 const auto loopSize = sumArray.extent(0);
943 using ExecSpaceType = typename DeviceType::execution_space;
944 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
945 Kokkos::parallel_for( policy, FunctorType(sumArray, inArray1, inArray2) );
946 }
947
948 // ------------------------------------------------------------------------------------
949
950 template<typename DeviceType>
951 template<typename inoutSumArrayValueType, class ...inoutSumArrayProperties,
952 typename inArrayValueType, class ...inArrayProperties>
953 void
954 RealSpaceTools<DeviceType>::
955 add( Kokkos::DynRankView<inoutSumArrayValueType,inoutSumArrayProperties...> inoutSumArray,
956 const Kokkos::DynRankView<inArrayValueType, inArrayProperties...> inArray ) {
957
958#ifdef HAVE_INTREPID2_DEBUG
959 {
960 INTREPID2_TEST_FOR_EXCEPTION( inArray.rank() != inoutSumArray.rank(), std::invalid_argument,
961 ">>> ERROR (RealSpaceTools::sum): Array arguments must have identical ranks!");
962 for (size_type i=0;i<inArray.rank();++i) {
963 INTREPID2_TEST_FOR_EXCEPTION( inArray.extent(i) != inoutSumArray.extent(i), std::invalid_argument,
964 ">>> ERROR (RealSpaceTools::sum): Dimensions of array arguments do not agree!");
965 }
966 }
967#endif
968 RealSpaceTools<DeviceType>::add(inoutSumArray, inoutSumArray, inArray);
969 }
970
971 // ------------------------------------------------------------------------------------
972
973 namespace FunctorRealSpaceTools {
977 template<typename diffArrayViewType,
978 typename inArray1ViewType,
979 typename inArray2ViewType>
980 struct F_subtract {
981 diffArrayViewType _diffArray;
982 const inArray1ViewType _inArray1;
983 const inArray2ViewType _inArray2;
984
985 KOKKOS_INLINE_FUNCTION
986 F_subtract( diffArrayViewType diffArray_,
987 inArray1ViewType inArray1_,
988 inArray2ViewType inArray2_ )
989 : _diffArray(diffArray_), _inArray1(inArray1_), _inArray2(inArray2_) {}
990
991 KOKKOS_INLINE_FUNCTION
992 void operator()(const ordinal_type i) const {
993 const ordinal_type jend = _diffArray.extent(1);
994 const ordinal_type kend = _diffArray.extent(2);
995 const ordinal_type lend = _diffArray.extent(3);
996 const ordinal_type mend = _diffArray.extent(4);
997
998 for (ordinal_type j=0;j<jend;++j)
999 for (ordinal_type k=0;k<kend;++k)
1000 for (ordinal_type l=0;l<lend;++l)
1001 for (ordinal_type m=0;m<mend;++m)
1002 _diffArray.access(i,j,k,l,m) = _inArray1.access(i,j,k,l,m) - _inArray2.access(i,j,k,l,m);
1003 }
1004 };
1005 }
1006
1007 template<typename DeviceType>
1008 template<typename diffArrayValueType, class ...diffArrayProperties,
1009 typename inArray1ValueType, class ...inArray1Properties,
1010 typename inArray2ValueType, class ...inArray2Properties>
1011 void
1012 RealSpaceTools<DeviceType>::
1013 subtract( Kokkos::DynRankView<diffArrayValueType,diffArrayProperties...> diffArray,
1014 const Kokkos::DynRankView<inArray1ValueType, inArray1Properties...> inArray1,
1015 const Kokkos::DynRankView<inArray2ValueType, inArray2Properties...> inArray2 ) {
1016
1017#ifdef HAVE_INTREPID2_DEBUG
1018 {
1019 INTREPID2_TEST_FOR_EXCEPTION( inArray1.rank() != inArray2.rank() ||
1020 inArray1.rank() != diffArray.rank(), std::invalid_argument,
1021 ">>> ERROR (RealSpaceTools::subtract): Array arguments must have identical ranks!");
1022 for (size_type i=0;i<inArray1.rank();++i) {
1023 INTREPID2_TEST_FOR_EXCEPTION( inArray1.extent(i) != inArray2.extent(i) ||
1024 inArray1.extent(i) != diffArray.extent(i), std::invalid_argument,
1025 ">>> ERROR (RealSpaceTools::subtract): Dimensions of array arguments do not agree!");
1026 }
1027 }
1028#endif
1029 using MemSpaceType = typename DeviceType::memory_space;
1030 constexpr bool are_accessible =
1031 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1032 typename decltype(diffArray)::memory_space>::accessible &&
1033 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1034 typename decltype(inArray1)::memory_space>::accessible &&
1035 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1036 typename decltype(inArray2)::memory_space>::accessible;
1037 static_assert(are_accessible, "RealSpaceTools<DeviceType>::subtract(..): input/output views' memory spaces are not compatible with DeviceType");
1038
1039 using FunctorType = FunctorRealSpaceTools::F_subtract<decltype(diffArray),decltype(inArray1),decltype(inArray2)>;
1040
1041 const size_type loopSize = diffArray.extent(0);
1042 using ExecSpaceType = typename DeviceType::execution_space;
1043 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
1044 Kokkos::parallel_for( policy, FunctorType(diffArray, inArray1, inArray2) );
1045 }
1046
1047 // ------------------------------------------------------------------------------------
1048
1049 template<typename DeviceType>
1050 template<typename inoutDiffArrayValueType, class ...inoutDiffArrayProperties,
1051 typename inArrayValueType, class ...inArrayProperties>
1052 void
1053 RealSpaceTools<DeviceType>::
1054 subtract( Kokkos::DynRankView<inoutDiffArrayValueType,inoutDiffArrayProperties...> inoutDiffArray,
1055 const Kokkos::DynRankView<inArrayValueType, inArrayProperties...> inArray ) {
1056#ifdef HAVE_INTREPID2_DEBUG
1057 {
1058 INTREPID2_TEST_FOR_EXCEPTION( inArray.rank() != inoutDiffArray.rank(), std::invalid_argument,
1059 ">>> ERROR (RealSpaceTools::subtract): Array arguments must have identical ranks!");
1060 for (size_type i=0;i<inArray.rank();++i) {
1061 INTREPID2_TEST_FOR_EXCEPTION( inArray.extent(i) != inoutDiffArray.extent(i), std::invalid_argument,
1062 ">>> ERROR (RealSpaceTools::subtract): Dimensions of array arguments do not agree!");
1063 }
1064 }
1065#endif
1066 RealSpaceTools<DeviceType>::subtract(inoutDiffArray, inoutDiffArray, inArray);
1067 }
1068
1069 // ------------------------------------------------------------------------------------
1070
1071 namespace FunctorRealSpaceTools {
1075 template<typename ValueType,
1076 typename scaledArrayViewType,
1077 typename inArrayViewType>
1078 struct F_scale {
1079 scaledArrayViewType _scaledArray;
1080 const inArrayViewType _inArray;
1081 const ValueType _alpha;
1082
1083 KOKKOS_INLINE_FUNCTION
1084 F_scale( scaledArrayViewType scaledArray_,
1085 inArrayViewType inArray_,
1086 const ValueType alpha_ )
1087 : _scaledArray(scaledArray_), _inArray(inArray_), _alpha(alpha_) {}
1088
1089 KOKKOS_INLINE_FUNCTION
1090 void operator()(const ordinal_type i) const {
1091 const ordinal_type jend = _inArray.extent(1);
1092 const ordinal_type kend = _inArray.extent(2);
1093 const ordinal_type lend = _inArray.extent(3);
1094 const ordinal_type mend = _inArray.extent(4);
1095
1096 for (ordinal_type j=0;j<jend;++j)
1097 for (ordinal_type k=0;k<kend;++k)
1098 for (ordinal_type l=0;l<lend;++l)
1099 for (ordinal_type m=0;m<mend;++m)
1100 _scaledArray.access(i,j,k,l,m) = _alpha*_inArray.access(i,j,k,l,m);
1101 }
1102 };
1103 }
1104
1105
1106 template<typename DeviceType>
1107 template<typename ValueType,
1108 typename scaledArrayValueType, class ...scaledArrayProperties,
1109 typename inArrayValueType, class ...inArrayProperties>
1110 void
1111 RealSpaceTools<DeviceType>::
1112 scale( Kokkos::DynRankView<scaledArrayValueType,scaledArrayProperties...> scaledArray,
1113 const Kokkos::DynRankView<inArrayValueType, inArrayProperties...> inArray,
1114 const ValueType alpha ) {
1115
1116#ifdef HAVE_INTREPID2_DEBUG
1117 {
1118 INTREPID2_TEST_FOR_EXCEPTION( inArray.rank() > 5, std::invalid_argument,
1119 ">>> ERROR (RealSpaceTools::scale): Input array container has rank larger than 5.");
1120 INTREPID2_TEST_FOR_EXCEPTION( inArray.rank() != scaledArray.rank(), std::invalid_argument,
1121 ">>> ERROR (RealSpaceTools::scale): Array arguments must have identical ranks!");
1122 for (size_type i=0;i<inArray.rank();++i) {
1123 INTREPID2_TEST_FOR_EXCEPTION( inArray.extent(i) != scaledArray.extent(i), std::invalid_argument,
1124 ">>> ERROR (RealSpaceTools::scale): Dimensions of array arguments do not agree!");
1125 }
1126 }
1127#endif
1128
1129 using MemSpaceType = typename DeviceType::memory_space;
1130 constexpr bool are_accessible =
1131 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1132 typename decltype(scaledArray)::memory_space>::accessible &&
1133 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1134 typename decltype(inArray)::memory_space>::accessible;
1135 static_assert(are_accessible, "RealSpaceTools<DeviceType>::scale(..): input/output views' memory spaces are not compatible with DeviceType");
1136
1137 using FunctorType = FunctorRealSpaceTools::F_scale<ValueType,decltype(scaledArray),decltype(inArray)>;
1138
1139 const auto loopSize = scaledArray.extent(0);
1140 using ExecSpaceType = typename DeviceType::execution_space;
1141 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
1142 Kokkos::parallel_for( policy, FunctorType(scaledArray, inArray, alpha) );
1143 }
1144
1145 // ------------------------------------------------------------------------------------
1146
1147 template<typename DeviceType>
1148 template<typename ValueType,
1149 typename inoutScaledArrayValueType, class ...inoutScaledArrayProperties>
1150 void
1151 RealSpaceTools<DeviceType>::
1152 scale( Kokkos::DynRankView<inoutScaledArrayValueType,inoutScaledArrayProperties...> inoutScaledArray,
1153 const ValueType alpha ) {
1154 RealSpaceTools<DeviceType>::scale(inoutScaledArray, inoutScaledArray, alpha);
1155 }
1156
1157
1158 // ------------------------------------------------------------------------------------
1159
1160 namespace FunctorRealSpaceTools {
1164 template<typename dotArrayViewType,
1165 typename inVec1ViewType,
1166 typename inVec2ViewType>
1167 struct F_dot {
1168 dotArrayViewType _dotArray;
1169 const inVec1ViewType _inVecs1;
1170 const inVec2ViewType _inVecs2;
1171
1172 KOKKOS_INLINE_FUNCTION
1173 F_dot( dotArrayViewType dotArray_,
1174 inVec1ViewType inVecs1_,
1175 inVec2ViewType inVecs2_ )
1176 : _dotArray(dotArray_), _inVecs1(inVecs1_), _inVecs2(inVecs2_) {}
1177
1178 KOKKOS_INLINE_FUNCTION
1179 void operator()(const ordinal_type iter) const {
1180 // the rank of normArray is either 1 or 2
1181 ordinal_type i, j;
1182 unrollIndex( i, j,
1183 _dotArray.extent(0),
1184 _dotArray.extent(1),
1185 iter );
1186
1187 const auto r = _inVecs1.rank();
1188 auto vec1 = ( r == 2 ? Kokkos::subview(_inVecs1, i, Kokkos::ALL()) :
1189 Kokkos::subview(_inVecs1, i, j, Kokkos::ALL()) );
1190 auto vec2 = ( r == 2 ? Kokkos::subview(_inVecs2, i, Kokkos::ALL()) :
1191 Kokkos::subview(_inVecs2, i, j, Kokkos::ALL()) );
1192
1193 _dotArray.access(i,j) = RealSpaceTools<>::Serial::dot(vec1, vec2);
1194 }
1195 };
1196 }
1197
1198 template<typename DeviceType>
1199 template<typename dotArrayValueType, class ...dotArrayProperties,
1200 typename inVec1ValueType, class ...inVec1Properties,
1201 typename inVec2ValueType, class ...inVec2Properties>
1202 void
1203 RealSpaceTools<DeviceType>::
1204 dot( Kokkos::DynRankView<dotArrayValueType,dotArrayProperties...> dotArray,
1205 const Kokkos::DynRankView<inVec1ValueType, inVec1Properties...> inVecs1,
1206 const Kokkos::DynRankView<inVec2ValueType, inVec2Properties...> inVecs2 ) {
1207
1208#ifdef HAVE_INTREPID2_DEBUG
1209 {
1210 INTREPID2_TEST_FOR_EXCEPTION( inVecs1.rank() != (dotArray.rank()+1), std::invalid_argument,
1211 ">>> ERROR (RealSpaceTools::dot): Ranks of norm and vector array arguments are incompatible!");
1212 INTREPID2_TEST_FOR_EXCEPTION( inVecs1.rank() != inVecs2.rank(), std::invalid_argument,
1213 ">>> ERROR (RealSpaceTools::dot): Ranks of input vector arguments must be identical!");
1214 INTREPID2_TEST_FOR_EXCEPTION( inVecs1.rank() < 2 || inVecs1.rank() > 3, std::invalid_argument,
1215 ">>> ERROR (RealSpaceTools::dot): Rank of input vector arguments must be 2 or 3!");
1216 for (size_type i=0;i<inVecs1.rank();++i) {
1217 INTREPID2_TEST_FOR_EXCEPTION( inVecs1.extent(i) != inVecs2.extent(i), std::invalid_argument,
1218 ">>> ERROR (RealSpaceTools::dot): Dimensions of input vector arguments do not agree!");
1219 }
1220 for (size_type i=0;i<dotArray.rank();++i) {
1221 INTREPID2_TEST_FOR_EXCEPTION( inVecs1.extent(i) != dotArray.extent(i), std::invalid_argument,
1222 ">>> ERROR (RealSpaceTools::dot): Dimensions of dot-product and vector arrays do not agree!");
1223 }
1224 }
1225#endif
1226 using MemSpaceType = typename DeviceType::memory_space;
1227 constexpr bool are_accessible =
1228 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1229 typename decltype(dotArray)::memory_space>::accessible &&
1230 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1231 typename decltype(inVecs1)::memory_space>::accessible &&
1232 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1233 typename decltype(inVecs2)::memory_space>::accessible;
1234 static_assert(are_accessible, "RealSpaceTools<DeviceType>::dot(..): input/output views' memory spaces are not compatible with DeviceType");
1235
1236 using FunctorType = FunctorRealSpaceTools::F_dot<decltype(dotArray),decltype(inVecs1),decltype(inVecs2)>;
1237
1238 const auto loopSize = dotArray.extent(0)*dotArray.extent(1);
1239 using ExecSpaceType = typename DeviceType::execution_space;
1240 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
1241 Kokkos::parallel_for( policy, FunctorType(dotArray, inVecs1, inVecs2) );
1242 }
1243
1244 // ------------------------------------------------------------------------------------
1245
1246 namespace FunctorRealSpaceTools {
1250 template<typename matVecViewType,
1251 typename inMatViewType,
1252 typename inVecViewType>
1253 struct F_matvec {
1254 matVecViewType _matVecs;
1255 const inMatViewType _inMats;
1256 const inVecViewType _inVecs;
1257
1258 KOKKOS_INLINE_FUNCTION
1259 F_matvec( matVecViewType matVecs_,
1260 inMatViewType inMats_,
1261 inVecViewType inVecs_ )
1262 : _matVecs(matVecs_), _inMats(inMats_), _inVecs(inVecs_) {}
1263
1264 KOKKOS_INLINE_FUNCTION
1265 void operator()(const ordinal_type iter) const {
1266 const ordinal_type rm = _inMats.rank(), rv = _inVecs.rank(), rr = _matVecs.rank();
1267 ordinal_type _i = iter, _j = 0;
1268
1269 if ( rr > 2 )
1270 unrollIndex( _i, _j,
1271 _matVecs.extent(0),
1272 _matVecs.extent(1),
1273 iter );
1274
1275 auto mat = ( rm == 2 ? Kokkos::subview(_inMats, Kokkos::ALL(), Kokkos::ALL()) :
1276 rm == 3 ? Kokkos::subview(_inMats, _i, Kokkos::ALL(), Kokkos::ALL()) :
1277 Kokkos::subview(_inMats, _i, _j, Kokkos::ALL(), Kokkos::ALL()) );
1278
1279 auto vec = ( rv == 1 ? Kokkos::subview(_inVecs, Kokkos::ALL()) :
1280 rv == 2 ? Kokkos::subview(_inVecs, _i, Kokkos::ALL()) :
1281 Kokkos::subview(_inVecs, _i, _j, Kokkos::ALL()) );
1282
1283 auto result = ( rr == 1 ? Kokkos::subview(_matVecs, Kokkos::ALL()) :
1284 rr == 2 ? Kokkos::subview(_matVecs, _i, Kokkos::ALL()) :
1285 Kokkos::subview(_matVecs, _i, _j, Kokkos::ALL()) );
1286
1287 const ordinal_type iend = result.extent(0);
1288 const ordinal_type jend = vec.extent(0);
1289
1290 for (ordinal_type i=0;i<iend;++i) {
1291 result(i) = 0;
1292 for (ordinal_type j=0;j<jend;++j)
1293 result(i) += mat(i, j)*vec(j);
1294 }
1295 }
1296 };
1297
1298
1302 template<typename outMatViewType,
1303 typename inMatViewType>
1304 struct F_AtA {
1305 outMatViewType _outMats;
1306 const inMatViewType _inMats;
1307
1308 KOKKOS_INLINE_FUNCTION
1309 F_AtA( outMatViewType outMats_,
1310 inMatViewType inMats_)
1311 : _outMats(outMats_), _inMats(inMats_) {}
1312
1313 KOKKOS_INLINE_FUNCTION
1314 void operator()(const ordinal_type iter) const {
1315 const ordinal_type rIM = _inMats.rank(), rOM = _outMats.rank();
1316 ordinal_type _i = iter, _j = 0;
1317
1318 if ( rOM > 3 )
1319 unrollIndex( _i, _j,
1320 _outMats.extent(0),
1321 _outMats.extent(1),
1322 iter );
1323
1324 auto inMat = ( rIM == 2 ? Kokkos::subview(_inMats, Kokkos::ALL(), Kokkos::ALL()) :
1325 rIM == 3 ? Kokkos::subview(_inMats, _i, Kokkos::ALL(), Kokkos::ALL()) :
1326 Kokkos::subview(_inMats, _i, _j, Kokkos::ALL(), Kokkos::ALL()) );
1327
1328 auto outMat = ( rOM == 2 ? Kokkos::subview(_outMats, Kokkos::ALL(), Kokkos::ALL()) :
1329 rOM == 3 ? Kokkos::subview(_outMats, _i, Kokkos::ALL(), Kokkos::ALL()) :
1330 Kokkos::subview(_outMats, _i, _j, Kokkos::ALL(), Kokkos::ALL()) );
1331
1332 const ordinal_type iend = outMat.extent(0);
1333 const ordinal_type jend = outMat.extent(1);
1334 const ordinal_type kend = inMat.extent(0);
1335
1336 for (ordinal_type i=0;i<iend;++i) {
1337 for (ordinal_type j=i;j<jend;++j) {
1338 outMat(i,j) = 0;
1339 for (ordinal_type k=0;k<kend;++k)
1340 outMat(i,j) += inMat(k, i)*inMat(k,j);
1341 }
1342 for (ordinal_type j=0;j<i;++j)
1343 outMat(i,j) = outMat(j,i);
1344 }
1345 }
1346 };
1347
1348 }
1349
1350 template<typename DeviceType>
1351 template<typename matVecValueType, class ...matVecProperties,
1352 typename inMatValueType, class ...inMatProperties,
1353 typename inVecValueType, class ...inVecProperties>
1354 void
1355 RealSpaceTools<DeviceType>::
1356 matvec( Kokkos::DynRankView<matVecValueType,matVecProperties...> matVecs,
1357 const Kokkos::DynRankView<inMatValueType, inMatProperties...> inMats,
1358 const Kokkos::DynRankView<inVecValueType, inVecProperties...> inVecs ) {
1359
1360#ifdef HAVE_INTREPID2_DEBUG
1361 INTREPID2_TEST_FOR_EXCEPTION( inMats.rank() < 2 || inMats.rank() > 4, std::invalid_argument,
1362 ">>> ERROR (RealSpaceTools::matvec): Rank of matrix array must be 2, 3 or 4!");
1363 INTREPID2_TEST_FOR_EXCEPTION( matVecs.rank() < 1 || matVecs.rank() > 3, std::invalid_argument,
1364 ">>> ERROR (RealSpaceTools::matvec): Rank of matVecs array must be 1, 2 or 3!");
1365 INTREPID2_TEST_FOR_EXCEPTION( inVecs.rank() < 1 || inVecs.rank() > 3, std::invalid_argument,
1366 ">>> ERROR (RealSpaceTools::matvec): Rank of inVecs array must be 1, 2 or 3!");
1367 if (inMats.rank() == 2) {
1368 // a single matrix, multiple input and output
1369 INTREPID2_TEST_FOR_EXCEPTION( matVecs.rank() != inVecs.rank(), std::invalid_argument,
1370 ">>> ERROR (RealSpaceTools::matvec): Vector arrays must be have the same rank!");
1371 // output must match
1372 for (ordinal_type i=0;i< (static_cast<ordinal_type>(inVecs.rank())-1);++i) {
1373 INTREPID2_TEST_FOR_EXCEPTION( matVecs.extent(i) != inVecs.extent(i), std::invalid_argument,
1374 ">>> ERROR (RealSpaceTools::matvec): Dimensions of vector array arguments do not agree!");
1375 }
1376 // matvec compatibility
1377 INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(0) != matVecs.extent(matVecs.rank()-1) ||
1378 inMats.extent(1) != inVecs.extent(inVecs.rank()-1), std::invalid_argument,
1379 ">>> ERROR (RealSpaceTools::matvec): Matvec dimensions are not compatible each other.");
1380 } else if (inVecs.rank() < matVecs.rank()) {
1381 // multiple matrix, single input and multiple output
1382 INTREPID2_TEST_FOR_EXCEPTION( inMats.rank() != (matVecs.rank()+1), std::invalid_argument,
1383 ">>> ERROR (RealSpaceTools::matvec): The result vector and matrix array arguments do not have compatible ranks!");
1384 for (ordinal_type i=0;i<(static_cast<ordinal_type>(inMats.rank())-2);++i) {
1385 INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(i) != matVecs.extent(i), std::invalid_argument,
1386 ">>> ERROR (RealSpaceTools::matvec): Dimensions of vector and matrix array arguments do not agree!");
1387 }
1388 // matvec compatibility
1389 INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(inMats.rank()-2) != matVecs.extent(matVecs.rank()-1) ||
1390 inMats.extent(inMats.rank()-1) != inVecs.extent(inVecs.rank()-1), std::invalid_argument,
1391 ">>> ERROR (RealSpaceTools::matvec): Matvec dimensions are not compatible each other.");
1392 } else {
1393 // multiple matrix, multiple input and multiple output
1394 INTREPID2_TEST_FOR_EXCEPTION( inMats.rank() != (matVecs.rank()+1), std::invalid_argument,
1395 ">>> ERROR (RealSpaceTools::matvec): The result vector and matrix array arguments do not have compatible ranks!");
1396 INTREPID2_TEST_FOR_EXCEPTION( matVecs.rank() != inVecs.rank(), std::invalid_argument,
1397 ">>> ERROR (RealSpaceTools::matvec): Vector arrays must be have the same rank!");
1398 for (ordinal_type i=0;i<(static_cast<ordinal_type>(inMats.rank())-2);++i) {
1399 INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(i) != matVecs.extent(i), std::invalid_argument,
1400 ">>> ERROR (RealSpaceTools::matvec): Dimensions of vector and matrix array arguments do not agree!");
1401 }
1402 for (size_type i=0;i<inVecs.rank();++i) {
1403 INTREPID2_TEST_FOR_EXCEPTION( matVecs.extent(i) != inVecs.extent(i), std::invalid_argument,
1404 ">>> ERROR (RealSpaceTools::matvec): Dimensions of vector array arguments do not agree!");
1405 }
1406 INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(inMats.rank()-1) != inVecs.extent(inVecs.rank()-1), std::invalid_argument,
1407 ">>> ERROR (RealSpaceTools::matvec): Matrix column dimension does not match to the length of a vector!");
1408 }
1409#endif
1410 using MemSpaceType = typename DeviceType::memory_space;
1411 constexpr bool are_accessible =
1412 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1413 typename decltype(matVecs)::memory_space>::accessible &&
1414 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1415 typename decltype(inMats)::memory_space>::accessible &&
1416 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1417 typename decltype(inVecs)::memory_space>::accessible;
1418 static_assert(are_accessible, "RealSpaceTools<DeviceType>::matvec(..): input/output views' memory spaces are not compatible with DeviceType");
1419
1420 using FunctorType = FunctorRealSpaceTools::F_matvec<decltype(matVecs),decltype(inMats),decltype(inVecs)>;
1421
1422 size_type loopSize = 1;
1423 const ordinal_type r = matVecs.rank() - 1;
1424 for (ordinal_type i=0;i<r;++i)
1425 loopSize *= matVecs.extent(i);
1426
1427 using ExecSpaceType = typename DeviceType::execution_space;
1428 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
1429 Kokkos::parallel_for( policy, FunctorType(matVecs, inMats, inVecs) );
1430 }
1431
1432
1433 template<typename DeviceType>
1434 template<typename outMatValueType, class ...outMatProperties,
1435 typename inMatValueType, class ...inMatProperties>
1436 void
1437 RealSpaceTools<DeviceType>::
1438 AtA( Kokkos::DynRankView<outMatValueType,outMatProperties...> outMats,
1439 const Kokkos::DynRankView<inMatValueType, inMatProperties...> inMats) {
1440
1441#ifdef HAVE_INTREPID2_DEBUG
1442 INTREPID2_TEST_FOR_EXCEPTION( inMats.rank() < 2 || inMats.rank() > 4, std::invalid_argument,
1443 ">>> ERROR (RealSpaceTools::AtA): Rank of input matrix array must be 2, 3 or 4!");
1444
1445 INTREPID2_TEST_FOR_EXCEPTION( inMats.rank() != outMats.rank(), std::invalid_argument,
1446 ">>> ERROR (RealSpaceTools::AtA): The matrices do not have the same ranks!");
1447 for (ordinal_type i=0;i<(static_cast<ordinal_type>(inMats.rank())-2);++i) {
1448 INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(i) != outMats.extent(i), std::invalid_argument,
1449 ">>> ERROR (RealSpaceTools::AtA): Dimensions of matrices arrays do not agree!");
1450 }
1451 // matmat compatibility
1452 INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(inMats.rank()-1) != outMats.extent(outMats.rank()-1) ||
1453 outMats.extent(outMats.rank()-1) != outMats.extent(outMats.rank()-2), std::invalid_argument,
1454 ">>> ERROR (RealSpaceTools::AtA): Matrices dimensions are not compatible.");
1455
1456#endif
1457 using MemSpaceType = typename DeviceType::memory_space;
1458 constexpr bool are_accessible =
1459 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1460 typename decltype(outMats)::memory_space>::accessible &&
1461 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1462 typename decltype(inMats)::memory_space>::accessible;
1463 static_assert(are_accessible, "RealSpaceTools<DeviceType>::AtA(..): input/output views' memory spaces are not compatible with DeviceType");
1464
1465 using FunctorType = FunctorRealSpaceTools::F_AtA<decltype(outMats),decltype(inMats)>;
1466
1467 size_type loopSize = 1;
1468 const ordinal_type r = outMats.rank() - 2;
1469 for (ordinal_type i=0;i<r;++i)
1470 loopSize *= outMats.extent(i);
1471
1472 using ExecSpaceType = typename DeviceType::execution_space;
1473 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
1474 Kokkos::parallel_for( policy, FunctorType(outMats, inMats) );
1475 }
1476
1477
1478 // ------------------------------------------------------------------------------------
1479
1480 namespace FunctorRealSpaceTools {
1484 template<typename vecProdViewType,
1485 typename inLeftViewType,
1486 typename inRightViewType>
1487 struct F_vecprod {
1488 vecProdViewType _vecProd;
1489 const inLeftViewType _inLeft;
1490 const inRightViewType _inRight;
1491 const bool _is_vecprod_3d;
1492
1493 KOKKOS_INLINE_FUNCTION
1494 F_vecprod( vecProdViewType vecProd_,
1495 inLeftViewType inLeft_,
1496 inRightViewType inRight_,
1497 const bool is_vecprod_3d_ )
1498 : _vecProd(vecProd_), _inLeft(inLeft_), _inRight(inRight_), _is_vecprod_3d(is_vecprod_3d_) {}
1499
1500 KOKKOS_INLINE_FUNCTION
1501 void operator()(const size_type iter) const {
1502 ordinal_type i, j;
1503 unrollIndex( i, j,
1504 _inLeft.extent(0),
1505 _inLeft.extent(1),
1506 iter );
1507
1508 const auto r = _inLeft.rank();
1509
1510 // branch prediction is possible
1511 if(r == 1) {
1512 if (_is_vecprod_3d) {
1513 _vecProd(0) = _inLeft(1)*_inRight(2) - _inLeft(2)*_inRight(1);
1514 _vecProd(1) = _inLeft(2)*_inRight(0) - _inLeft(0)*_inRight(2);
1515 _vecProd(2) = _inLeft(0)*_inRight(1) - _inLeft(1)*_inRight(0);
1516 } else {
1517 _vecProd(0) = _inLeft(0)*_inRight(1) - _inLeft(1)*_inRight(0);
1518 }
1519 } else if(r == 2) {
1520 if (_is_vecprod_3d) {
1521 _vecProd(i, 0) = _inLeft(i, 1)*_inRight(i, 2) - _inLeft(i, 2)*_inRight(i, 1);
1522 _vecProd(i, 1) = _inLeft(i, 2)*_inRight(i, 0) - _inLeft(i, 0)*_inRight(i, 2);
1523 _vecProd(i, 2) = _inLeft(i, 0)*_inRight(i, 1) - _inLeft(i, 1)*_inRight(i, 0);
1524 } else {
1525 _vecProd(i, 0) = _inLeft(i, 0)*_inRight(i, 1) - _inLeft(i, 1)*_inRight(i, 0);
1526 }
1527 } else {
1528 if (_is_vecprod_3d) {
1529 _vecProd(i, j, 0) = _inLeft(i, j, 1)*_inRight(i, j, 2) - _inLeft(i, j, 2)*_inRight(i, j, 1);
1530 _vecProd(i, j, 1) = _inLeft(i, j, 2)*_inRight(i, j, 0) - _inLeft(i, j, 0)*_inRight(i, j, 2);
1531 _vecProd(i, j, 2) = _inLeft(i, j, 0)*_inRight(i, j, 1) - _inLeft(i, j, 1)*_inRight(i, j, 0);
1532 } else {
1533 _vecProd(i, j, 0) = _inLeft(i, j, 0)*_inRight(i, j, 1) - _inLeft(i, j, 1)*_inRight(i, j, 0);
1534 }
1535 }
1536 }
1537 };
1538 }
1539
1540 template<typename DeviceType>
1541 template<typename vecProdValueType, class ...vecProdProperties,
1542 typename inLeftValueType, class ...inLeftProperties,
1543 typename inRightValueType, class ...inRightProperties>
1544 void
1545 RealSpaceTools<DeviceType>::
1546 vecprod( Kokkos::DynRankView<vecProdValueType,vecProdProperties...> vecProd,
1547 const Kokkos::DynRankView<inLeftValueType, inLeftProperties...> inLeft,
1548 const Kokkos::DynRankView<inRightValueType,inRightProperties...> inRight ) {
1549
1550#ifdef HAVE_INTREPID2_DEBUG
1551 {
1552 /*
1553 * Check array rank and spatial dimension range (if applicable)
1554 * (1) all array arguments are required to have matching dimensions and rank: (D), (I0,D) or (I0,I1,D)
1555 * (2) spatial dimension should be 2 or 3
1556 */
1557 // (1) check rank range on inLeft and then compare the other arrays with inLeft
1558 INTREPID2_TEST_FOR_EXCEPTION( inLeft.rank() < 1 || inLeft.rank() > 3, std::invalid_argument,
1559 ">>> ERROR (RealSpaceTools::vecprod): Rank of matrix array must be 1, 2, or 3!");
1560 INTREPID2_TEST_FOR_EXCEPTION( inLeft.rank() != inRight.rank(), std::invalid_argument,
1561 ">>> ERROR (RealSpaceTools::vecprod): Right and left arrays must be have the same rank!");
1562 INTREPID2_TEST_FOR_EXCEPTION( inLeft.rank() != vecProd.rank(), std::invalid_argument,
1563 ">>> ERROR (RealSpaceTools::vecprod): Left and vecProd arrays must be have the same rank!");
1564 for (size_type i=0;i<inLeft.rank();++i) {
1565 INTREPID2_TEST_FOR_EXCEPTION( inLeft.extent(i) != inRight.extent(i), std::invalid_argument,
1566 ">>> ERROR (RealSpaceTools::vecprod): Dimensions of matrix arguments do not agree!");
1567 INTREPID2_TEST_FOR_EXCEPTION( inLeft.extent(i) != vecProd.extent(i), std::invalid_argument,
1568 ">>> ERROR (RealSpaceTools::vecprod): Dimensions of matrix arguments do not agree!");
1569 }
1570
1571 // (2) spatial dimension ordinal = array rank - 1. Suffices to check only one array because we just
1572 // checked whether or not the arrays have matching dimensions.
1573 INTREPID2_TEST_FOR_EXCEPTION( inLeft.extent(inLeft.rank()-1) < 2 ||
1574 inLeft.extent(inLeft.rank()-1) > 3, std::invalid_argument,
1575 ">>> ERROR (RealSpaceTools::vecprod): Dimensions of arrays (rank-1) must be 2 or 3!");
1576 }
1577#endif
1578 using MemSpaceType = typename DeviceType::memory_space;
1579 constexpr bool are_accessible =
1580 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1581 typename decltype(vecProd)::memory_space>::accessible &&
1582 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1583 typename decltype(inLeft)::memory_space>::accessible &&
1584 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1585 typename decltype(inRight)::memory_space>::accessible;
1586 static_assert(are_accessible, "RealSpaceTools<DeviceType>::vecprod(..): input/output views' memory spaces are not compatible with DeviceType");
1587
1588 using FunctorType = FunctorRealSpaceTools::F_vecprod<decltype(vecProd),decltype(inLeft),decltype(inRight)>;
1589
1590 const auto r = inLeft.rank();
1591 const auto loopSize = ( r == 1 ? 1 :
1592 r == 2 ? inLeft.extent(0) :
1593 inLeft.extent(0)*inLeft.extent(1) );
1594 const bool is_vecprod_3d = (inLeft.extent(inLeft.rank() - 1) == 3);
1595 using ExecSpaceType = typename DeviceType::execution_space;
1596 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
1597 Kokkos::parallel_for( policy, FunctorType(vecProd, inLeft, inRight, is_vecprod_3d) );
1598 }
1599
1600 // ------------------------------------------------------------------------------------
1601
1602} // namespace Intrepid2
1603
1604#endif
1605
1606
1607
1608
1609
1610// =====================================
1611// Too much things...
1612//
1613
1614// template<class ...inVec1Properties,
1615// class ...inVec2Properties>
1616// KOKKOS_INLINE_FUNCTION
1617// typename Kokkos::DynRankView<inVec1Properties...>::value_type
1618// RealSpaceTools<DeviceType>::
1619// dot( const Kokkos::DynRankView<inVec1Properties...> inVec1,
1620// const Kokkos::DynRankView<inVec2Properties...> inVec2 ) {
1621
1622// #ifdef HAVE_INTREPID2_DEBUG
1623// INTREPID2_TEST_FOR_ABORT( inVec1.rank != 1 || inVec2.rank() != 1,
1624// ">>> ERROR (RealSpaceTools::dot): Vector arguments must have rank 1!");
1625// INTREPID2_TEST_FOR_ABORT( inVec1.extent(0) != inVec2.extent(0),
1626// ">>> ERROR (RealSpaceTools::dot): Dimensions of vector arguments must agree!");
1627// #endif
1628// typedef typename Kokkos::DynRankView<inVec1Properties...>::value_type value_type;
1629
1630// // designed for small problems
1631// value_type r_val(0);
1632
1633// // * This is probably not necessary
1634// if (std::is_same<ExecSpace,Kokkos::Serial>::value) {
1635// const ordinal_type iend = iVec1.extent(0);
1636// for (ordinal_type i=0;i<iend;++i)
1637// r_val += inVec1(i)*inVec2(i);
1638// } else {
1639// struct Functor {
1640// Kokkos::DynRankView<inVec1Properties...> _inVec1;
1641// Kokkos::DynRankView<inVec2Properties...> _inVec2;
1642
1643// KOKKOS_INLINE_FUNCTION
1644// Functor(const Kokkos::DynRankView<inVec1Properties...> inVec1_,
1645// const Kokkos::DynRankView<inVec2Properties...> inVec2_);
1646
1647// KOKKOS_INLINE_FUNCTION
1648// void operator()(ordinal_type i, value_type &dst ) const {
1649// dst += _inVec1(i)*_inVec2(i);
1650// }
1651
1652// KOKKOS_INLINE_FUNCTION
1653// void join(volatile value_type &dst ,
1654// const volatile value_type &src) const {
1655// dst += src;
1656// }
1657// };
1658// const ordinal_type iend = inVec1.extent(0);
1659// Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, iend);
1660// Kokkos::parallel_for( policy, Functor(inVec1, inVec2), r_val );
1661// }
1662// return r_val;
1663// }
1664
1665
Intrepid2::getFunctorRank
enable_if_t< has_rank_method< Functor >::value, unsigned > KOKKOS_INLINE_FUNCTION getFunctorRank(const Functor &functor)
Definition Intrepid2_Utils.hpp:747
Intrepid2::get_scalar_value
KOKKOS_FORCEINLINE_FUNCTION constexpr std::enable_if<!(std::is_standard_layout< T >::value &&std::is_trivial< T >::value), typenameScalarTraits< T >::scalar_type >::type get_scalar_value(const T &obj)
functions returning the scalar value. for pod types, they return the input object itself....
Definition Intrepid2_Utils.hpp:286
Intrepid2::RealSpaceTools
Implementation of basic linear algebra functionality in Euclidean space.
Definition Intrepid2_RealSpaceTools.hpp:46
Intrepid2::RealSpaceTools::extractScalarValues
static void extractScalarValues(Kokkos::DynRankView< outputValueType, outputProperties... > output, const Kokkos::DynRankView< inputValueType, inputProperties... > input)
Extract scalar type values from Sacado-based array.
Definition Intrepid2_RealSpaceToolsDef.hpp:224
Intrepid2::RealSpaceTools::add
static void add(Kokkos::DynRankView< sumArrayValueType, sumArrayProperties... > sumArray, const Kokkos::DynRankView< inArray1ValueType, inArray1Properties... > inArray1, const Kokkos::DynRankView< inArray2ValueType, inArray2Properties... > inArray2)
Adds arrays inArray1 and inArray2: sumArray = inArray1 + inArray2.
Definition Intrepid2_RealSpaceToolsDef.hpp:914
Intrepid2::RealSpaceTools::vectorNorm
static void vectorNorm(Kokkos::DynRankView< normArrayValueType, normArrayProperties... > normArray, const Kokkos::DynRankView< inVecValueType, inVecProperties... > inVecs, const ENorm normType)
Computes norms (1, 2, infinity) of vectors stored in a array of total rank 2 (array of vectors),...
Definition Intrepid2_RealSpaceToolsDef.hpp:467
Intrepid2::RealSpaceTools::inverse
static void inverse(InverseMatrixViewType inverseMats, MatrixViewType inMats)
Computes inverses of nonsingular matrices stored in an array of total rank 2 (single matrix),...
Definition Intrepid2_RealSpaceToolsDef.hpp:717
Intrepid2::RealSpaceTools::matvec
static void matvec(Kokkos::DynRankView< matVecValueType, matVecProperties... > matVecs, const Kokkos::DynRankView< inMatValueType, inMatProperties... > inMats, const Kokkos::DynRankView< inVecValueType, inVecProperties... > inVecs)
Matrix-vector left multiply using multidimensional arrays: matVec = inMat * inVec.
Definition Intrepid2_RealSpaceToolsDef.hpp:1356
Intrepid2::RealSpaceTools::absval
static void absval(Kokkos::DynRankView< absArrayValueType, absArrayProperties... > absArray, const Kokkos::DynRankView< inArrayValueType, inArrayProperties... > inArray)
Computes absolute value of an array.
Definition Intrepid2_RealSpaceToolsDef.hpp:386
Intrepid2::RealSpaceTools::AtA
static void AtA(Kokkos::DynRankView< outMatValueType, outMatProperties... > outMats, const Kokkos::DynRankView< inMatValueType, inMatProperties... > inMats)
Computes the matrix-matrix product , for an input rectangular matrix .
Definition Intrepid2_RealSpaceToolsDef.hpp:1438
Intrepid2::RealSpaceTools::vecprod
static void vecprod(Kokkos::DynRankView< vecProdValueType, vecProdProperties... > vecProd, const Kokkos::DynRankView< inLeftValueType, inLeftProperties... > inLeft, const Kokkos::DynRankView< inRightValueType, inRightProperties... > inRight)
Vector product using multidimensional arrays: vecProd = inVecLeft x inVecRight
Definition Intrepid2_RealSpaceToolsDef.hpp:1546
Intrepid2::RealSpaceTools::subtract
static void subtract(Kokkos::DynRankView< diffArrayValueType, diffArrayProperties... > diffArray, const Kokkos::DynRankView< inArray1ValueType, inArray1Properties... > inArray1, const Kokkos::DynRankView< inArray2ValueType, inArray2Properties... > inArray2)
Subtracts inArray2 from inArray1: diffArray = inArray1 - inArray2.
Definition Intrepid2_RealSpaceToolsDef.hpp:1013
Intrepid2::RealSpaceTools::scale
static void scale(Kokkos::DynRankView< scaledArrayValueType, scaledArrayProperties... > scaledArray, const Kokkos::DynRankView< inArrayValueType, inArrayProperties... > inArray, const ValueType alpha)
Multiplies array inArray by the scalar scalar (componentwise): scaledArray = scalar * inArray.
Definition Intrepid2_RealSpaceToolsDef.hpp:1112
Intrepid2::RealSpaceTools::dot
static void dot(Kokkos::DynRankView< dotArrayValueType, dotArrayProperties... > dotArray, const Kokkos::DynRankView< inVec1ValueType, inVec1Properties... > inVecs1, const Kokkos::DynRankView< inVec2ValueType, inVec2Properties... > inVecs2)
Computes dot product of vectors stored in an array of total rank 2 (array of vectors),...
Definition Intrepid2_RealSpaceToolsDef.hpp:1204
Intrepid2::RealSpaceTools::det
static void det(DeterminantArrayViewType detArray, const MatrixViewType inMats)
Computes determinants of matrices stored in an array of total rank 3 (array of matrices),...
Definition Intrepid2_RealSpaceToolsDef.hpp:823
Intrepid2::RealSpaceTools::transpose
static void transpose(Kokkos::DynRankView< transposeMatValueType, transposeMatProperties... > transposeMats, const Kokkos::DynRankView< inMatValueType, inMatProperties... > inMats)
Computes transposes of square matrices stored in an array of total rank 2 (single matrix),...
Definition Intrepid2_RealSpaceToolsDef.hpp:552
Intrepid2::RealSpaceTools::clone
static void clone(Kokkos::DynRankView< outputValueType, outputProperties... > output, const Kokkos::DynRankView< inputValueType, inputProperties... > input)
Clone input array.
Definition Intrepid2_RealSpaceToolsDef.hpp:301
Intrepid2::Util
small utility functions
Definition Intrepid2_Utils.hpp:238
Intrepid2::FunctorRealSpaceTools::F_AtA
Functor to compute matvec see Intrepid2::RealSpaceTools for more.
Definition Intrepid2_RealSpaceToolsDef.hpp:1304
Intrepid2::FunctorRealSpaceTools::F_absval
Functor to compute absolute value see Intrepid2::RealSpaceTools for more.
Definition Intrepid2_RealSpaceToolsDef.hpp:354
Intrepid2::FunctorRealSpaceTools::F_add
Functor to add md arrays see Intrepid2::RealSpaceTools for more.
Definition Intrepid2_RealSpaceToolsDef.hpp:881
Intrepid2::FunctorRealSpaceTools::F_clone
Functor for clone see Intrepid2::RealSpaceTools for more.
Definition Intrepid2_RealSpaceToolsDef.hpp:249
Intrepid2::FunctorRealSpaceTools::F_det
Functor to compute determinant see Intrepid2::RealSpaceTools for more.
Definition Intrepid2_RealSpaceToolsDef.hpp:769
Intrepid2::FunctorRealSpaceTools::F_dot
Functor to compute dot product see Intrepid2::RealSpaceTools for more.
Definition Intrepid2_RealSpaceToolsDef.hpp:1167
Intrepid2::FunctorRealSpaceTools::F_extractScalarValues
Functor for extractScalarValues see Intrepid2::RealSpaceTools for more.
Definition Intrepid2_RealSpaceToolsDef.hpp:194
Intrepid2::FunctorRealSpaceTools::F_inverse
Functor to compute inverse see Intrepid2::RealSpaceTools for more.
Definition Intrepid2_RealSpaceToolsDef.hpp:598
Intrepid2::FunctorRealSpaceTools::F_matvec
Functor to compute matvec see Intrepid2::RealSpaceTools for more.
Definition Intrepid2_RealSpaceToolsDef.hpp:1253
Intrepid2::FunctorRealSpaceTools::F_scale
Functor to scale md arrays see Intrepid2::RealSpaceTools for more.
Definition Intrepid2_RealSpaceToolsDef.hpp:1078
Intrepid2::FunctorRealSpaceTools::F_subtract
Functor to subtract md arrays see Intrepid2::RealSpaceTools for more.
Definition Intrepid2_RealSpaceToolsDef.hpp:980
Intrepid2::FunctorRealSpaceTools::F_transpose
Functor to compute transpose see Intrepid2::RealSpaceTools for more.
Definition Intrepid2_RealSpaceToolsDef.hpp:507
Intrepid2::FunctorRealSpaceTools::F_vecprod
Functor to compute vecprod see Intrepid2::RealSpaceTools for more.
Definition Intrepid2_RealSpaceToolsDef.hpp:1487
Intrepid2::FunctorRealSpaceTools::F_vectorNorm
Functor to compute vector norm see Intrepid2::RealSpaceTools for more.
Definition Intrepid2_RealSpaceToolsDef.hpp:435
Intrepid2::RealSpaceTools::Serial::det
static KOKKOS_INLINE_FUNCTION MatrixViewType::value_type det(const MatrixViewType inMat)
Computes determinant of a single square matrix stored in an array of rank 2.
Definition Intrepid2_RealSpaceToolsDef.hpp:104
Intrepid2::RealSpaceTools::Serial::vectorNorm
static KOKKOS_INLINE_FUNCTION inVecValueType vectorNorm(const Kokkos::DynRankView< inVecValueType, inVecProperties... > inVec, const ENorm normType)
Computes norm (1, 2, infinity) of a single vector stored in an array of rank 1.
Definition Intrepid2_RealSpaceToolsDef.hpp:34
Intrepid2::RealSpaceTools::Serial::dot
static KOKKOS_INLINE_FUNCTION inVec1ValueType dot(const Kokkos::DynRankView< inVec1ValueType, inVec1Properties... > inVec1, const Kokkos::DynRankView< inVec2ValueType, inVec2Properties... > inVec2)
Computes dot product of two vectors stored in arrays of rank 1.
Definition Intrepid2_RealSpaceToolsDef.hpp:152
Generated by doxygen 1.9.8