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