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