Intrepid2
Intrepid2_RealSpaceToolsDef.hpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Intrepid Package
5// Copyright (2007) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Kyungjoo Kim (kyukim@sandia.gov), or
38// Mauro Perego (mperego@sandia.gov)
39//
40// ************************************************************************
41// @HEADER
42
49#ifndef __INTREPID2_REALSPACETOOLS_DEF_HPP__
50#define __INTREPID2_REALSPACETOOLS_DEF_HPP__
51
52namespace Intrepid2 {
53
54 // ------------------------------------------------------------------------------------
55
56 //
57 // serial version code
58 //
59
60 // ------------------------------------------------------------------------------------
61
62 template<typename DeviceType>
63 template<typename inVecValueType, class ...inVecProperties>
64 KOKKOS_INLINE_FUNCTION
65 inVecValueType
67 vectorNorm( const Kokkos::DynRankView<inVecValueType,inVecProperties...> inVec,
68 const ENorm normType ) {
69#ifdef HAVE_INTREPID2_DEBUG
70 {
71 bool dbgInfo = false;
72 INTREPID2_TEST_FOR_DEBUG_ABORT( !(inVec.rank() >= 1 && inVec.rank() <= 5), dbgInfo,
73 ">>> ERROR (RealSpaceTools::vectorNorm): Vector argument must have rank 1!" );
74#ifdef INTREPID2_TEST_FOR_DEBUG_ABORT_OVERRIDE_TO_CONTINUE
75 if (dbgInfo) return inVecValueType(0);
76#endif
77 }
78#endif
79 typedef inVecValueType value_type;
80
81 value_type norm(0);
82
83 const ordinal_type iend = inVec.extent(0);
84 const ordinal_type jend = inVec.extent(1);
85 const ordinal_type kend = inVec.extent(2);
86 const ordinal_type lend = inVec.extent(3);
87 const ordinal_type mend = inVec.extent(4);
88 switch(normType) {
89 case NORM_TWO:{
90 for (ordinal_type i=0;i<iend;++i)
91 for (ordinal_type j=0;j<jend;++j)
92 for (ordinal_type k=0;k<kend;++k)
93 for (ordinal_type l=0;l<lend;++l)
94 for (ordinal_type m=0;m<mend;++m)
95 norm += inVec.access(i,j,k,l,m)*inVec.access(i,j,k,l,m);
96 norm = sqrt(norm);
97 break;
98 }
99 case NORM_INF:{
100 for (ordinal_type i=0;i<iend;++i)
101 for (ordinal_type j=0;j<jend;++j)
102 for (ordinal_type k=0;k<kend;++k)
103 for (ordinal_type l=0;l<lend;++l)
104 for (ordinal_type m=0;m<mend;++m) {
105 const value_type current = Util<value_type>::abs(inVec.access(i,j,k,l,m));
106 norm = (norm < current ? current : norm);
107 }
108 break;
109 }
110 case NORM_ONE:{
111 for (ordinal_type i=0;i<iend;++i)
112 for (ordinal_type j=0;j<jend;++j)
113 for (ordinal_type k=0;k<kend;++k)
114 for (ordinal_type l=0;l<lend;++l)
115 for (ordinal_type m=0;m<mend;++m)
116 norm += Util<value_type>::abs(inVec.access(i,j,k,l,m));
117 break;
118 }
119 default: {
120 INTREPID2_TEST_FOR_ABORT( ( (normType != NORM_TWO) &&
121 (normType != NORM_INF) &&
122 (normType != NORM_ONE) ),
123 ">>> ERROR (RealSpaceTools::vectorNorm): Invalid argument normType.");
124 }
125 }
126 return norm;
127 }
128
129 // ------------------------------------------------------------------------------------
130
131
132 template<typename DeviceType>
133 template<class MatrixViewType>
134 KOKKOS_INLINE_FUNCTION
135 typename MatrixViewType::value_type
137 det( const MatrixViewType inMat ) {
138
139 typedef typename decltype(inMat)::non_const_value_type value_type;
140#ifdef HAVE_INTREPID2_DEBUG
141 {
142 bool dbgInfo = false;
143 INTREPID2_TEST_FOR_DEBUG_ABORT( getFunctorRank( inMat ) != 2, dbgInfo,
144 ">>> ERROR (RealSpaceTools::det): Rank of matrix argument must be 2!");
145 INTREPID2_TEST_FOR_DEBUG_ABORT( inMat.extent(0) != inMat.extent(1), dbgInfo,
146 ">>> ERROR (RealSpaceTools::det): Matrix is not square!");
147 INTREPID2_TEST_FOR_DEBUG_ABORT( inMat.extent(0) < 1 || inMat.extent(0) > 3, dbgInfo,
148 ">>> ERROR (RealSpaceTools::det): Spatial dimension must be 1, 2, or 3!");
149#ifdef INTREPID2_TEST_FOR_DEBUG_ABORT_OVERRIDE_TO_CONTINUE
150 if (dbgInfo) return value_type(0);
151#endif
152 }
153#endif
154 const auto dim = inMat.extent(0);
155
156 value_type r_val = 0.0;
157 switch (dim) {
158 case 3:
159 r_val = ( inMat(0,0) * inMat(1,1) * inMat(2,2) +
160 inMat(1,0) * inMat(2,1) * inMat(0,2) +
161 inMat(2,0) * inMat(0,1) * inMat(1,2) -
162 inMat(2,0) * inMat(1,1) * inMat(0,2) -
163 inMat(0,0) * inMat(2,1) * inMat(1,2) -
164 inMat(1,0) * inMat(0,1) * inMat(2,2) );
165 break;
166 case 2:
167 r_val = ( inMat(0,0) * inMat(1,1) -
168 inMat(0,1) * inMat(1,0) );
169 break;
170 case 1:
171 r_val = ( inMat(0,0) );
172 break;
173 }
174 return r_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
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>
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(1);
239 const ordinal_type kend = _input.extent(2);
240 const ordinal_type lend = _input.extent(3);
241 const ordinal_type mend = _input.extent(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
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
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
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
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>
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(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
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
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<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(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 switch (mat.extent(0)) {
663 case 1: {
664 inv(0,0) = value_type(1)/mat(0,0);
665 break;
666 }
667 case 2: {
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 break;
674 }
675 case 3: {
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 break;
702 }
703 }
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<class InverseMatrixViewType, class MatrixViewType>
748 void
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;
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 // ------------------------------------------------------------------------------------
795
796 namespace FunctorRealSpaceTools {
800 template<typename detArrayViewType,
801 typename inMatViewType, int rank>
802 struct F_det {
803 detArrayViewType _detArray;
804 inMatViewType _inMats;
805
806 KOKKOS_INLINE_FUNCTION
807 F_det( detArrayViewType detArray_,
808 inMatViewType inMats_ )
809 : _detArray(detArray_), _inMats(inMats_) {}
810
811 template< bool B, class T = void >
812 using enable_if_t = typename std::enable_if<B,T>::type;
813
814 template<int M=rank>
815 KOKKOS_INLINE_FUNCTION
816 enable_if_t<M==1 && supports_rank_3<inMatViewType>::value >
817 operator()(const ordinal_type pt) const {
818 const auto mat = Kokkos::subview(_inMats, pt, Kokkos::ALL(), Kokkos::ALL());
819 _detArray(pt) = RealSpaceTools<>::Serial::det(mat);
820 }
821
822 template<int M=rank>
823 KOKKOS_INLINE_FUNCTION
824 enable_if_t<M==1 && !supports_rank_3<inMatViewType>::value >
825 operator()(const ordinal_type pt) const {
826 // empty implementation to allow compilation of parallel_for. (this combination never gets invoked at run-time, but does occur at compile-time.)
827 }
828
829 template<int M=rank>
830 KOKKOS_INLINE_FUNCTION
831 enable_if_t<M==2 && supports_rank_4<inMatViewType>::value >
832 operator()(const ordinal_type cl,
833 const ordinal_type pt) const {
834 const auto mat = Kokkos::subview(_inMats, cl, pt, Kokkos::ALL(), Kokkos::ALL());
835 _detArray(cl, pt) = RealSpaceTools<>::Serial::det(mat);
836 }
837
838 template<int M=rank>
839 KOKKOS_INLINE_FUNCTION
840 enable_if_t<M==2 && !supports_rank_4<inMatViewType>::value >
841 operator()(const ordinal_type cl,
842 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 }
847
848template<class DeterminantArrayViewType, class MatrixViewType>
849static void
850det( DeterminantArrayViewType detArray, const MatrixViewType inMats );
851
852 template<typename DeviceType>
853 template<class DeterminantArrayViewType, class MatrixViewType>
854 void
856 det( DeterminantArrayViewType detArray, const MatrixViewType inMats ) {
857
858#ifdef HAVE_INTREPID2_DEBUG
859 {
860 const unsigned matrixRank = getFunctorRank(inMats);
861 const unsigned detRank = getFunctorRank(detArray);
862 INTREPID2_TEST_FOR_EXCEPTION( matrixRank != detRank+2, std::invalid_argument,
863 ">>> ERROR (RealSpaceTools::det): Determinant and matrix array arguments do not have compatible ranks!");
864 INTREPID2_TEST_FOR_EXCEPTION( matrixRank < 3 || matrixRank > 4, std::invalid_argument,
865 ">>> ERROR (RealSpaceTools::det): Rank of matrix array must be 3 or 4!");
866 for (size_type i=0;i<matrixRank-2;++i) {
867 INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(i) != detArray.extent(i), std::invalid_argument,
868 ">>> ERROR (RealSpaceTools::det): Dimensions of determinant and matrix array arguments do not agree!");
869 }
870 INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(matrixRank-2) != inMats.extent(matrixRank-1), std::invalid_argument,
871 ">>> ERROR (RealSpaceTools::det): Matrices are not square!");
872 INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(matrixRank-2) < 1 || inMats.extent(matrixRank-2) > 3, std::invalid_argument,
873 ">>> ERROR (RealSpaceTools::det): Spatial dimension must be 1, 2, or 3!");
874 }
875#endif
876
877 typedef typename DeviceType::execution_space ExecSpaceType;
878
879 const int detArrayRank = getFunctorRank(detArray);
880
881 switch (detArrayRank) {
882 case 1: { // output P and input P,D,D
884 using range_policy_type = Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> >;
885 range_policy_type policy(0, detArray.extent(0));
886 Kokkos::parallel_for( policy, FunctorType(detArray, inMats) );
887 break;
888 }
889 case 2: { // output C,P and input C,P,D,D
891 using range_policy_type = Kokkos::MDRangePolicy
892 < ExecSpaceType, Kokkos::Rank<2>, Kokkos::IndexType<ordinal_type> >;
893 range_policy_type policy( { 0, 0 },
894 { detArray.extent(0), detArray.extent(1) } );
895 Kokkos::parallel_for( policy, FunctorType(detArray, inMats) );
896 break;
897 }
898 default: {
899 INTREPID2_TEST_FOR_EXCEPTION( true, std::invalid_argument,
900 ">>> ERROR (RealSpaceTools::det): Rank of detArray must be 1 or 2!");
901 }
902 }
903 }
904
905 // ------------------------------------------------------------------------------------
906
907 namespace FunctorRealSpaceTools {
911 template<typename sumArrayViewType,
912 typename inArray1Viewtype,
913 typename inArray2ViewType>
914 struct F_add {
915 sumArrayViewType _sumArray;
916 inArray1Viewtype _inArray1;
917 inArray2ViewType _inArray2;
918
919 KOKKOS_INLINE_FUNCTION
920 F_add( sumArrayViewType sumArray_,
921 inArray1Viewtype inArray1_,
922 inArray2ViewType inArray2_ )
923 : _sumArray(sumArray_), _inArray1(inArray1_), _inArray2(inArray2_) {}
924
925 KOKKOS_INLINE_FUNCTION
926 void operator()(const ordinal_type i) const {
927 const ordinal_type jend = _sumArray.extent(1);
928 const ordinal_type kend = _sumArray.extent(2);
929 const ordinal_type lend = _sumArray.extent(3);
930 const ordinal_type mend = _sumArray.extent(4);
931
932 for (ordinal_type j=0;j<jend;++j)
933 for (ordinal_type k=0;k<kend;++k)
934 for (ordinal_type l=0;l<lend;++l)
935 for (ordinal_type m=0;m<mend;++m)
936 _sumArray.access(i,j,k,l,m) = _inArray1.access(i,j,k,l,m) + _inArray2.access(i,j,k,l,m);
937 }
938 };
939 }
940
941 template<typename DeviceType>
942 template<typename sumArrayValueType, class ...sumArrayProperties,
943 typename inArray1ValueType, class ...inArray1Properties,
944 typename inArray2ValueType, class ...inArray2Properties>
945 void
947 add( Kokkos::DynRankView<sumArrayValueType,sumArrayProperties...> sumArray,
948 const Kokkos::DynRankView<inArray1ValueType,inArray1Properties...> inArray1,
949 const Kokkos::DynRankView<inArray2ValueType,inArray2Properties...> inArray2 ) {
950
951#ifdef HAVE_INTREPID2_DEBUG
952 {
953 INTREPID2_TEST_FOR_EXCEPTION( inArray1.rank() != inArray2.rank() ||
954 inArray1.rank() != sumArray.rank(), std::invalid_argument,
955 ">>> ERROR (RealSpaceTools::add): Array arguments must have identical ranks!");
956 for (size_type i=0;i<inArray1.rank();++i) {
957 INTREPID2_TEST_FOR_EXCEPTION( inArray1.extent(i) != inArray2.extent(i) ||
958 inArray1.extent(i) != sumArray.extent(i), std::invalid_argument,
959 ">>> ERROR (RealSpaceTools::add): Dimensions of array arguments do not agree!");
960 }
961 }
962#endif
963 using MemSpaceType = typename DeviceType::memory_space;
964 constexpr bool are_accessible =
965 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
966 typename decltype(sumArray)::memory_space>::accessible &&
967 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
968 typename decltype(inArray1)::memory_space>::accessible &&
969 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
970 typename decltype(inArray2)::memory_space>::accessible;
971 static_assert(are_accessible, "RealSpaceTools<DeviceType>::add(..): input/output views' memory spaces are not compatible with DeviceType");
972
973 using FunctorType = FunctorRealSpaceTools::F_add<decltype(sumArray),decltype(inArray1),decltype(inArray2)>;
974
975 const auto loopSize = sumArray.extent(0);
976 using ExecSpaceType = typename DeviceType::execution_space;
977 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
978 Kokkos::parallel_for( policy, FunctorType(sumArray, inArray1, inArray2) );
979 }
980
981 // ------------------------------------------------------------------------------------
982
983 template<typename DeviceType>
984 template<typename inoutSumArrayValueType, class ...inoutSumArrayProperties,
985 typename inArrayValueType, class ...inArrayProperties>
986 void
988 add( Kokkos::DynRankView<inoutSumArrayValueType,inoutSumArrayProperties...> inoutSumArray,
989 const Kokkos::DynRankView<inArrayValueType, inArrayProperties...> inArray ) {
990
991#ifdef HAVE_INTREPID2_DEBUG
992 {
993 INTREPID2_TEST_FOR_EXCEPTION( inArray.rank() != inoutSumArray.rank(), std::invalid_argument,
994 ">>> ERROR (RealSpaceTools::sum): Array arguments must have identical ranks!");
995 for (size_type i=0;i<inArray.rank();++i) {
996 INTREPID2_TEST_FOR_EXCEPTION( inArray.extent(i) != inoutSumArray.extent(i), std::invalid_argument,
997 ">>> ERROR (RealSpaceTools::sum): Dimensions of array arguments do not agree!");
998 }
999 }
1000#endif
1001 RealSpaceTools<DeviceType>::add(inoutSumArray, inoutSumArray, inArray);
1002 }
1003
1004 // ------------------------------------------------------------------------------------
1005
1006 namespace FunctorRealSpaceTools {
1010 template<typename diffArrayViewType,
1011 typename inArray1ViewType,
1012 typename inArray2ViewType>
1013 struct F_subtract {
1014 diffArrayViewType _diffArray;
1015 const inArray1ViewType _inArray1;
1016 const inArray2ViewType _inArray2;
1017
1018 KOKKOS_INLINE_FUNCTION
1019 F_subtract( diffArrayViewType diffArray_,
1020 inArray1ViewType inArray1_,
1021 inArray2ViewType inArray2_ )
1022 : _diffArray(diffArray_), _inArray1(inArray1_), _inArray2(inArray2_) {}
1023
1024 KOKKOS_INLINE_FUNCTION
1025 void operator()(const ordinal_type i) const {
1026 const ordinal_type jend = _diffArray.extent(1);
1027 const ordinal_type kend = _diffArray.extent(2);
1028 const ordinal_type lend = _diffArray.extent(3);
1029 const ordinal_type mend = _diffArray.extent(4);
1030
1031 for (ordinal_type j=0;j<jend;++j)
1032 for (ordinal_type k=0;k<kend;++k)
1033 for (ordinal_type l=0;l<lend;++l)
1034 for (ordinal_type m=0;m<mend;++m)
1035 _diffArray.access(i,j,k,l,m) = _inArray1.access(i,j,k,l,m) - _inArray2.access(i,j,k,l,m);
1036 }
1037 };
1038 }
1039
1040 template<typename DeviceType>
1041 template<typename diffArrayValueType, class ...diffArrayProperties,
1042 typename inArray1ValueType, class ...inArray1Properties,
1043 typename inArray2ValueType, class ...inArray2Properties>
1044 void
1046 subtract( Kokkos::DynRankView<diffArrayValueType,diffArrayProperties...> diffArray,
1047 const Kokkos::DynRankView<inArray1ValueType, inArray1Properties...> inArray1,
1048 const Kokkos::DynRankView<inArray2ValueType, inArray2Properties...> inArray2 ) {
1049
1050#ifdef HAVE_INTREPID2_DEBUG
1051 {
1052 INTREPID2_TEST_FOR_EXCEPTION( inArray1.rank() != inArray2.rank() ||
1053 inArray1.rank() != diffArray.rank(), std::invalid_argument,
1054 ">>> ERROR (RealSpaceTools::subtract): Array arguments must have identical ranks!");
1055 for (size_type i=0;i<inArray1.rank();++i) {
1056 INTREPID2_TEST_FOR_EXCEPTION( inArray1.extent(i) != inArray2.extent(i) ||
1057 inArray1.extent(i) != diffArray.extent(i), std::invalid_argument,
1058 ">>> ERROR (RealSpaceTools::subtract): Dimensions of array arguments do not agree!");
1059 }
1060 }
1061#endif
1062 using MemSpaceType = typename DeviceType::memory_space;
1063 constexpr bool are_accessible =
1064 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1065 typename decltype(diffArray)::memory_space>::accessible &&
1066 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1067 typename decltype(inArray1)::memory_space>::accessible &&
1068 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1069 typename decltype(inArray2)::memory_space>::accessible;
1070 static_assert(are_accessible, "RealSpaceTools<DeviceType>::subtract(..): input/output views' memory spaces are not compatible with DeviceType");
1071
1072 using FunctorType = FunctorRealSpaceTools::F_subtract<decltype(diffArray),decltype(inArray1),decltype(inArray2)>;
1073
1074 const size_type loopSize = diffArray.extent(0);
1075 using ExecSpaceType = typename DeviceType::execution_space;
1076 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
1077 Kokkos::parallel_for( policy, FunctorType(diffArray, inArray1, inArray2) );
1078 }
1079
1080 // ------------------------------------------------------------------------------------
1081
1082 template<typename DeviceType>
1083 template<typename inoutDiffArrayValueType, class ...inoutDiffArrayProperties,
1084 typename inArrayValueType, class ...inArrayProperties>
1085 void
1087 subtract( Kokkos::DynRankView<inoutDiffArrayValueType,inoutDiffArrayProperties...> inoutDiffArray,
1088 const Kokkos::DynRankView<inArrayValueType, inArrayProperties...> inArray ) {
1089#ifdef HAVE_INTREPID2_DEBUG
1090 {
1091 INTREPID2_TEST_FOR_EXCEPTION( inArray.rank() != inoutDiffArray.rank(), std::invalid_argument,
1092 ">>> ERROR (RealSpaceTools::subtract): Array arguments must have identical ranks!");
1093 for (size_type i=0;i<inArray.rank();++i) {
1094 INTREPID2_TEST_FOR_EXCEPTION( inArray.extent(i) != inoutDiffArray.extent(i), std::invalid_argument,
1095 ">>> ERROR (RealSpaceTools::subtract): Dimensions of array arguments do not agree!");
1096 }
1097 }
1098#endif
1099 RealSpaceTools<DeviceType>::subtract(inoutDiffArray, inoutDiffArray, inArray);
1100 }
1101
1102 // ------------------------------------------------------------------------------------
1103
1104 namespace FunctorRealSpaceTools {
1108 template<typename ValueType,
1109 typename scaledArrayViewType,
1110 typename inArrayViewType>
1111 struct F_scale {
1112 scaledArrayViewType _scaledArray;
1113 const inArrayViewType _inArray;
1114 const ValueType _alpha;
1115
1116 KOKKOS_INLINE_FUNCTION
1117 F_scale( scaledArrayViewType scaledArray_,
1118 inArrayViewType inArray_,
1119 const ValueType alpha_ )
1120 : _scaledArray(scaledArray_), _inArray(inArray_), _alpha(alpha_) {}
1121
1122 KOKKOS_INLINE_FUNCTION
1123 void operator()(const ordinal_type i) const {
1124 const ordinal_type jend = _inArray.extent(1);
1125 const ordinal_type kend = _inArray.extent(2);
1126 const ordinal_type lend = _inArray.extent(3);
1127 const ordinal_type mend = _inArray.extent(4);
1128
1129 for (ordinal_type j=0;j<jend;++j)
1130 for (ordinal_type k=0;k<kend;++k)
1131 for (ordinal_type l=0;l<lend;++l)
1132 for (ordinal_type m=0;m<mend;++m)
1133 _scaledArray.access(i,j,k,l,m) = _alpha*_inArray.access(i,j,k,l,m);
1134 }
1135 };
1136 }
1137
1138
1139 template<typename DeviceType>
1140 template<typename ValueType,
1141 typename scaledArrayValueType, class ...scaledArrayProperties,
1142 typename inArrayValueType, class ...inArrayProperties>
1143 void
1145 scale( Kokkos::DynRankView<scaledArrayValueType,scaledArrayProperties...> scaledArray,
1146 const Kokkos::DynRankView<inArrayValueType, inArrayProperties...> inArray,
1147 const ValueType alpha ) {
1148
1149#ifdef HAVE_INTREPID2_DEBUG
1150 {
1151 INTREPID2_TEST_FOR_EXCEPTION( inArray.rank() > 5, std::invalid_argument,
1152 ">>> ERROR (RealSpaceTools::scale): Input array container has rank larger than 5.");
1153 INTREPID2_TEST_FOR_EXCEPTION( inArray.rank() != scaledArray.rank(), std::invalid_argument,
1154 ">>> ERROR (RealSpaceTools::scale): Array arguments must have identical ranks!");
1155 for (size_type i=0;i<inArray.rank();++i) {
1156 INTREPID2_TEST_FOR_EXCEPTION( inArray.extent(i) != scaledArray.extent(i), std::invalid_argument,
1157 ">>> ERROR (RealSpaceTools::scale): Dimensions of array arguments do not agree!");
1158 }
1159 }
1160#endif
1161
1162 using MemSpaceType = typename DeviceType::memory_space;
1163 constexpr bool are_accessible =
1164 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1165 typename decltype(scaledArray)::memory_space>::accessible &&
1166 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1167 typename decltype(inArray)::memory_space>::accessible;
1168 static_assert(are_accessible, "RealSpaceTools<DeviceType>::scale(..): input/output views' memory spaces are not compatible with DeviceType");
1169
1170 using FunctorType = FunctorRealSpaceTools::F_scale<ValueType,decltype(scaledArray),decltype(inArray)>;
1171
1172 const auto loopSize = scaledArray.extent(0);
1173 using ExecSpaceType = typename DeviceType::execution_space;
1174 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
1175 Kokkos::parallel_for( policy, FunctorType(scaledArray, inArray, alpha) );
1176 }
1177
1178 // ------------------------------------------------------------------------------------
1179
1180 template<typename DeviceType>
1181 template<typename ValueType,
1182 typename inoutScaledArrayValueType, class ...inoutScaledArrayProperties>
1183 void
1185 scale( Kokkos::DynRankView<inoutScaledArrayValueType,inoutScaledArrayProperties...> inoutScaledArray,
1186 const ValueType alpha ) {
1187 RealSpaceTools<DeviceType>::scale(inoutScaledArray, inoutScaledArray, alpha);
1188 }
1189
1190
1191 // ------------------------------------------------------------------------------------
1192
1193 namespace FunctorRealSpaceTools {
1197 template<typename dotArrayViewType,
1198 typename inVec1ViewType,
1199 typename inVec2ViewType>
1200 struct F_dot {
1201 dotArrayViewType _dotArray;
1202 const inVec1ViewType _inVecs1;
1203 const inVec2ViewType _inVecs2;
1204
1205 KOKKOS_INLINE_FUNCTION
1206 F_dot( dotArrayViewType dotArray_,
1207 inVec1ViewType inVecs1_,
1208 inVec2ViewType inVecs2_ )
1209 : _dotArray(dotArray_), _inVecs1(inVecs1_), _inVecs2(inVecs2_) {}
1210
1211 KOKKOS_INLINE_FUNCTION
1212 void operator()(const ordinal_type iter) const {
1213 // the rank of normArray is either 1 or 2
1214 ordinal_type i, j;
1215 unrollIndex( i, j,
1216 _dotArray.extent(0),
1217 _dotArray.extent(1),
1218 iter );
1219
1220 const auto r = _inVecs1.rank();
1221 auto vec1 = ( r == 2 ? Kokkos::subview(_inVecs1, i, Kokkos::ALL()) :
1222 Kokkos::subview(_inVecs1, i, j, Kokkos::ALL()) );
1223 auto vec2 = ( r == 2 ? Kokkos::subview(_inVecs2, i, Kokkos::ALL()) :
1224 Kokkos::subview(_inVecs2, i, j, Kokkos::ALL()) );
1225
1226 _dotArray(i,j) = RealSpaceTools<>::Serial::dot(vec1, vec2);
1227 }
1228 };
1229 }
1230
1231 template<typename DeviceType>
1232 template<typename dotArrayValueType, class ...dotArrayProperties,
1233 typename inVec1ValueType, class ...inVec1Properties,
1234 typename inVec2ValueType, class ...inVec2Properties>
1235 void
1237 dot( Kokkos::DynRankView<dotArrayValueType,dotArrayProperties...> dotArray,
1238 const Kokkos::DynRankView<inVec1ValueType, inVec1Properties...> inVecs1,
1239 const Kokkos::DynRankView<inVec2ValueType, inVec2Properties...> inVecs2 ) {
1240
1241#ifdef HAVE_INTREPID2_DEBUG
1242 {
1243 INTREPID2_TEST_FOR_EXCEPTION( inVecs1.rank() != (dotArray.rank()+1), std::invalid_argument,
1244 ">>> ERROR (RealSpaceTools::dot): Ranks of norm and vector array arguments are incompatible!");
1245 INTREPID2_TEST_FOR_EXCEPTION( inVecs1.rank() != inVecs2.rank(), std::invalid_argument,
1246 ">>> ERROR (RealSpaceTools::dot): Ranks of input vector arguments must be identical!");
1247 INTREPID2_TEST_FOR_EXCEPTION( inVecs1.rank() < 2 || inVecs1.rank() > 3, std::invalid_argument,
1248 ">>> ERROR (RealSpaceTools::dot): Rank of input vector arguments must be 2 or 3!");
1249 for (size_type i=0;i<inVecs1.rank();++i) {
1250 INTREPID2_TEST_FOR_EXCEPTION( inVecs1.extent(i) != inVecs2.extent(i), std::invalid_argument,
1251 ">>> ERROR (RealSpaceTools::dot): Dimensions of input vector arguments do not agree!");
1252 }
1253 for (size_type i=0;i<dotArray.rank();++i) {
1254 INTREPID2_TEST_FOR_EXCEPTION( inVecs1.extent(i) != dotArray.extent(i), std::invalid_argument,
1255 ">>> ERROR (RealSpaceTools::dot): Dimensions of dot-product and vector arrays do not agree!");
1256 }
1257 }
1258#endif
1259 using MemSpaceType = typename DeviceType::memory_space;
1260 constexpr bool are_accessible =
1261 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1262 typename decltype(dotArray)::memory_space>::accessible &&
1263 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1264 typename decltype(inVecs1)::memory_space>::accessible &&
1265 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1266 typename decltype(inVecs2)::memory_space>::accessible;
1267 static_assert(are_accessible, "RealSpaceTools<DeviceType>::dot(..): input/output views' memory spaces are not compatible with DeviceType");
1268
1269 using FunctorType = FunctorRealSpaceTools::F_dot<decltype(dotArray),decltype(inVecs1),decltype(inVecs2)>;
1270
1271 const auto loopSize = dotArray.extent(0)*dotArray.extent(1);
1272 using ExecSpaceType = typename DeviceType::execution_space;
1273 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
1274 Kokkos::parallel_for( policy, FunctorType(dotArray, inVecs1, inVecs2) );
1275 }
1276
1277 // ------------------------------------------------------------------------------------
1278
1279 namespace FunctorRealSpaceTools {
1283 template<typename matVecViewType,
1284 typename inMatViewType,
1285 typename inVecViewType>
1286 struct F_matvec {
1287 matVecViewType _matVecs;
1288 const inMatViewType _inMats;
1289 const inVecViewType _inVecs;
1290
1291 KOKKOS_INLINE_FUNCTION
1292 F_matvec( matVecViewType matVecs_,
1293 inMatViewType inMats_,
1294 inVecViewType inVecs_ )
1295 : _matVecs(matVecs_), _inMats(inMats_), _inVecs(inVecs_) {}
1296
1297 KOKKOS_INLINE_FUNCTION
1298 void operator()(const ordinal_type iter) const {
1299 const ordinal_type rm = _inMats.rank(), rv = _inVecs.rank(), rr = _matVecs.rank();
1300 ordinal_type _i = iter, _j = 0;
1301
1302 if ( rr > 2 )
1303 unrollIndex( _i, _j,
1304 _matVecs.extent(0),
1305 _matVecs.extent(1),
1306 iter );
1307
1308 auto mat = ( rm == 2 ? Kokkos::subview(_inMats, Kokkos::ALL(), Kokkos::ALL()) :
1309 rm == 3 ? Kokkos::subview(_inMats, _i, Kokkos::ALL(), Kokkos::ALL()) :
1310 Kokkos::subview(_inMats, _i, _j, Kokkos::ALL(), Kokkos::ALL()) );
1311
1312 auto vec = ( rv == 1 ? Kokkos::subview(_inVecs, Kokkos::ALL()) :
1313 rv == 2 ? Kokkos::subview(_inVecs, _i, Kokkos::ALL()) :
1314 Kokkos::subview(_inVecs, _i, _j, Kokkos::ALL()) );
1315
1316 auto result = ( rr == 1 ? Kokkos::subview(_matVecs, Kokkos::ALL()) :
1317 rr == 2 ? Kokkos::subview(_matVecs, _i, Kokkos::ALL()) :
1318 Kokkos::subview(_matVecs, _i, _j, Kokkos::ALL()) );
1319
1320 const ordinal_type iend = result.extent(0);
1321 const ordinal_type jend = vec.extent(0);
1322
1323 for (ordinal_type i=0;i<iend;++i) {
1324 result(i) = 0;
1325 for (ordinal_type j=0;j<jend;++j)
1326 result(i) += mat(i, j)*vec(j);
1327 }
1328 }
1329 };
1330
1331
1335 template<typename outMatViewType,
1336 typename inMatViewType>
1337 struct F_AtA {
1338 outMatViewType _outMats;
1339 const inMatViewType _inMats;
1340
1341 KOKKOS_INLINE_FUNCTION
1342 F_AtA( outMatViewType outMats_,
1343 inMatViewType inMats_)
1344 : _outMats(outMats_), _inMats(inMats_) {}
1345
1346 KOKKOS_INLINE_FUNCTION
1347 void operator()(const ordinal_type iter) const {
1348 const ordinal_type rIM = _inMats.rank(), rOM = _outMats.rank();
1349 ordinal_type _i = iter, _j = 0;
1350
1351 if ( rOM > 3 )
1352 unrollIndex( _i, _j,
1353 _outMats.extent(0),
1354 _outMats.extent(1),
1355 iter );
1356
1357 auto inMat = ( rIM == 2 ? Kokkos::subview(_inMats, Kokkos::ALL(), Kokkos::ALL()) :
1358 rIM == 3 ? Kokkos::subview(_inMats, _i, Kokkos::ALL(), Kokkos::ALL()) :
1359 Kokkos::subview(_inMats, _i, _j, Kokkos::ALL(), Kokkos::ALL()) );
1360
1361 auto outMat = ( rOM == 2 ? Kokkos::subview(_outMats, Kokkos::ALL(), Kokkos::ALL()) :
1362 rOM == 3 ? Kokkos::subview(_outMats, _i, Kokkos::ALL(), Kokkos::ALL()) :
1363 Kokkos::subview(_outMats, _i, _j, Kokkos::ALL(), Kokkos::ALL()) );
1364
1365 const ordinal_type iend = outMat.extent(0);
1366 const ordinal_type jend = outMat.extent(1);
1367 const ordinal_type kend = inMat.extent(0);
1368
1369 for (ordinal_type i=0;i<iend;++i) {
1370 for (ordinal_type j=i;j<jend;++j) {
1371 outMat(i,j) = 0;
1372 for (ordinal_type k=0;k<kend;++k)
1373 outMat(i,j) += inMat(k, i)*inMat(k,j);
1374 }
1375 for (ordinal_type j=0;j<i;++j)
1376 outMat(i,j) = outMat(j,i);
1377 }
1378 }
1379 };
1380
1381 }
1382
1383 template<typename DeviceType>
1384 template<typename matVecValueType, class ...matVecProperties,
1385 typename inMatValueType, class ...inMatProperties,
1386 typename inVecValueType, class ...inVecProperties>
1387 void
1389 matvec( Kokkos::DynRankView<matVecValueType,matVecProperties...> matVecs,
1390 const Kokkos::DynRankView<inMatValueType, inMatProperties...> inMats,
1391 const Kokkos::DynRankView<inVecValueType, inVecProperties...> inVecs ) {
1392
1393#ifdef HAVE_INTREPID2_DEBUG
1394 INTREPID2_TEST_FOR_EXCEPTION( inMats.rank() < 2 || inMats.rank() > 4, std::invalid_argument,
1395 ">>> ERROR (RealSpaceTools::matvec): Rank of matrix array must be 2, 3 or 4!");
1396 INTREPID2_TEST_FOR_EXCEPTION( matVecs.rank() < 1 || matVecs.rank() > 3, std::invalid_argument,
1397 ">>> ERROR (RealSpaceTools::matvec): Rank of matVecs array must be 1, 2 or 3!");
1398 INTREPID2_TEST_FOR_EXCEPTION( inVecs.rank() < 1 || inVecs.rank() > 3, std::invalid_argument,
1399 ">>> ERROR (RealSpaceTools::matvec): Rank of inVecs array must be 1, 2 or 3!");
1400 if (inMats.rank() == 2) {
1401 // a single matrix, multiple input and output
1402 INTREPID2_TEST_FOR_EXCEPTION( matVecs.rank() != inVecs.rank(), std::invalid_argument,
1403 ">>> ERROR (RealSpaceTools::matvec): Vector arrays must be have the same rank!");
1404 // output must match
1405 for (ordinal_type i=0;i< (static_cast<ordinal_type>(inVecs.rank())-1);++i) {
1406 INTREPID2_TEST_FOR_EXCEPTION( matVecs.extent(i) != inVecs.extent(i), std::invalid_argument,
1407 ">>> ERROR (RealSpaceTools::matvec): Dimensions of vector array arguments do not agree!");
1408 }
1409 // matvec compatibility
1410 INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(0) != matVecs.extent(matVecs.rank()-1) ||
1411 inMats.extent(1) != inVecs.extent(inVecs.rank()-1), std::invalid_argument,
1412 ">>> ERROR (RealSpaceTools::matvec): Matvec dimensions are not compatible each other.");
1413 } else if (inVecs.rank() < matVecs.rank()) {
1414 // multiple matrix, single input and multiple output
1415 INTREPID2_TEST_FOR_EXCEPTION( inMats.rank() != (matVecs.rank()+1), std::invalid_argument,
1416 ">>> ERROR (RealSpaceTools::matvec): The result vector and matrix array arguments do not have compatible ranks!");
1417 for (ordinal_type i=0;i<(static_cast<ordinal_type>(inMats.rank())-2);++i) {
1418 INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(i) != matVecs.extent(i), std::invalid_argument,
1419 ">>> ERROR (RealSpaceTools::matvec): Dimensions of vector and matrix array arguments do not agree!");
1420 }
1421 // matvec compatibility
1422 INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(inMats.rank()-2) != matVecs.extent(matVecs.rank()-1) ||
1423 inMats.extent(inMats.rank()-1) != inVecs.extent(inVecs.rank()-1), std::invalid_argument,
1424 ">>> ERROR (RealSpaceTools::matvec): Matvec dimensions are not compatible each other.");
1425 } else {
1426 // multiple matrix, multiple input and multiple output
1427 INTREPID2_TEST_FOR_EXCEPTION( inMats.rank() != (matVecs.rank()+1), std::invalid_argument,
1428 ">>> ERROR (RealSpaceTools::matvec): The result vector and matrix array arguments do not have compatible ranks!");
1429 INTREPID2_TEST_FOR_EXCEPTION( matVecs.rank() != inVecs.rank(), std::invalid_argument,
1430 ">>> ERROR (RealSpaceTools::matvec): Vector arrays must be have the same rank!");
1431 for (ordinal_type i=0;i<(static_cast<ordinal_type>(inMats.rank())-2);++i) {
1432 INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(i) != matVecs.extent(i), std::invalid_argument,
1433 ">>> ERROR (RealSpaceTools::matvec): Dimensions of vector and matrix array arguments do not agree!");
1434 }
1435 for (size_type i=0;i<inVecs.rank();++i) {
1436 INTREPID2_TEST_FOR_EXCEPTION( matVecs.extent(i) != inVecs.extent(i), std::invalid_argument,
1437 ">>> ERROR (RealSpaceTools::matvec): Dimensions of vector array arguments do not agree!");
1438 }
1439 INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(inMats.rank()-1) != inVecs.extent(inVecs.rank()-1), std::invalid_argument,
1440 ">>> ERROR (RealSpaceTools::matvec): Matrix column dimension does not match to the length of a vector!");
1441 }
1442#endif
1443 using MemSpaceType = typename DeviceType::memory_space;
1444 constexpr bool are_accessible =
1445 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1446 typename decltype(matVecs)::memory_space>::accessible &&
1447 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1448 typename decltype(inMats)::memory_space>::accessible &&
1449 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1450 typename decltype(inVecs)::memory_space>::accessible;
1451 static_assert(are_accessible, "RealSpaceTools<DeviceType>::matvec(..): input/output views' memory spaces are not compatible with DeviceType");
1452
1453 using FunctorType = FunctorRealSpaceTools::F_matvec<decltype(matVecs),decltype(inMats),decltype(inVecs)>;
1454
1455 size_type loopSize = 1;
1456 const ordinal_type r = matVecs.rank() - 1;
1457 for (ordinal_type i=0;i<r;++i)
1458 loopSize *= matVecs.extent(i);
1459
1460 using ExecSpaceType = typename DeviceType::execution_space;
1461 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
1462 Kokkos::parallel_for( policy, FunctorType(matVecs, inMats, inVecs) );
1463 }
1464
1465
1466 template<typename DeviceType>
1467 template<typename outMatValueType, class ...outMatProperties,
1468 typename inMatValueType, class ...inMatProperties>
1469 void
1471 AtA( Kokkos::DynRankView<outMatValueType,outMatProperties...> outMats,
1472 const Kokkos::DynRankView<inMatValueType, inMatProperties...> inMats) {
1473
1474#ifdef HAVE_INTREPID2_DEBUG
1475 INTREPID2_TEST_FOR_EXCEPTION( inMats.rank() < 2 || inMats.rank() > 4, std::invalid_argument,
1476 ">>> ERROR (RealSpaceTools::AtA): Rank of input matrix array must be 2, 3 or 4!");
1477
1478 INTREPID2_TEST_FOR_EXCEPTION( inMats.rank() != outMats.rank(), std::invalid_argument,
1479 ">>> ERROR (RealSpaceTools::AtA): The matrices do not have the same ranks!");
1480 for (ordinal_type i=0;i<(static_cast<ordinal_type>(inMats.rank())-2);++i) {
1481 INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(i) != outMats.extent(i), std::invalid_argument,
1482 ">>> ERROR (RealSpaceTools::AtA): Dimensions of matrices arrays do not agree!");
1483 }
1484 // matmat compatibility
1485 INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(inMats.rank()-1) != outMats.extent(outMats.rank()-1) ||
1486 outMats.extent(outMats.rank()-1) != outMats.extent(outMats.rank()-2), std::invalid_argument,
1487 ">>> ERROR (RealSpaceTools::AtA): Matrices dimensions are not compatible.");
1488
1489#endif
1490 using MemSpaceType = typename DeviceType::memory_space;
1491 constexpr bool are_accessible =
1492 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1493 typename decltype(outMats)::memory_space>::accessible &&
1494 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1495 typename decltype(inMats)::memory_space>::accessible;
1496 static_assert(are_accessible, "RealSpaceTools<DeviceType>::AtA(..): input/output views' memory spaces are not compatible with DeviceType");
1497
1498 using FunctorType = FunctorRealSpaceTools::F_AtA<decltype(outMats),decltype(inMats)>;
1499
1500 size_type loopSize = 1;
1501 const ordinal_type r = outMats.rank() - 2;
1502 for (ordinal_type i=0;i<r;++i)
1503 loopSize *= outMats.extent(i);
1504
1505 using ExecSpaceType = typename DeviceType::execution_space;
1506 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
1507 Kokkos::parallel_for( policy, FunctorType(outMats, inMats) );
1508 }
1509
1510
1511 // ------------------------------------------------------------------------------------
1512
1513 namespace FunctorRealSpaceTools {
1517 template<typename vecProdViewType,
1518 typename inLeftViewType,
1519 typename inRightViewType>
1520 struct F_vecprod {
1521 vecProdViewType _vecProd;
1522 const inLeftViewType _inLeft;
1523 const inRightViewType _inRight;
1524 const bool _is_vecprod_3d;
1525
1526 KOKKOS_INLINE_FUNCTION
1527 F_vecprod( vecProdViewType vecProd_,
1528 inLeftViewType inLeft_,
1529 inRightViewType inRight_,
1530 const bool is_vecprod_3d_ )
1531 : _vecProd(vecProd_), _inLeft(inLeft_), _inRight(inRight_), _is_vecprod_3d(is_vecprod_3d_) {}
1532
1533 KOKKOS_INLINE_FUNCTION
1534 void operator()(const size_type iter) const {
1535 ordinal_type i, j;
1536 unrollIndex( i, j,
1537 _inLeft.extent(0),
1538 _inLeft.extent(1),
1539 iter );
1540
1541 const auto r = _inLeft.rank();
1542
1543 auto left = ( r == 1 ? Kokkos::subview(_inLeft, Kokkos::ALL()) :
1544 r == 2 ? Kokkos::subview(_inLeft, i, Kokkos::ALL()) :
1545 Kokkos::subview(_inLeft, i, j, Kokkos::ALL()) );
1546
1547 auto right = ( r == 1 ? Kokkos::subview(_inRight, Kokkos::ALL()) :
1548 r == 2 ? Kokkos::subview(_inRight, i, Kokkos::ALL()) :
1549 Kokkos::subview(_inRight, i, j, Kokkos::ALL()) );
1550
1551 auto result = ( r == 1 ? Kokkos::subview(_vecProd, Kokkos::ALL()) :
1552 r == 2 ? Kokkos::subview(_vecProd, i, Kokkos::ALL()) :
1553 Kokkos::subview(_vecProd, i, j, Kokkos::ALL()) );
1554
1555 // branch prediction is possible
1556 if (_is_vecprod_3d) {
1557 result(0) = left(1)*right(2) - left(2)*right(1);
1558 result(1) = left(2)*right(0) - left(0)*right(2);
1559 result(2) = left(0)*right(1) - left(1)*right(0);
1560 } else {
1561 result(0) = left(0)*right(1) - left(1)*right(0);
1562 }
1563 }
1564 };
1565 }
1566
1567 template<typename DeviceType>
1568 template<typename vecProdValueType, class ...vecProdProperties,
1569 typename inLeftValueType, class ...inLeftProperties,
1570 typename inRightValueType, class ...inRightProperties>
1571 void
1573 vecprod( Kokkos::DynRankView<vecProdValueType,vecProdProperties...> vecProd,
1574 const Kokkos::DynRankView<inLeftValueType, inLeftProperties...> inLeft,
1575 const Kokkos::DynRankView<inRightValueType,inRightProperties...> inRight ) {
1576
1577#ifdef HAVE_INTREPID2_DEBUG
1578 {
1579 /*
1580 * Check array rank and spatial dimension range (if applicable)
1581 * (1) all array arguments are required to have matching dimensions and rank: (D), (I0,D) or (I0,I1,D)
1582 * (2) spatial dimension should be 2 or 3
1583 */
1584 // (1) check rank range on inLeft and then compare the other arrays with inLeft
1585 INTREPID2_TEST_FOR_EXCEPTION( inLeft.rank() < 1 || inLeft.rank() > 3, std::invalid_argument,
1586 ">>> ERROR (RealSpaceTools::vecprod): Rank of matrix array must be 1, 2, or 3!");
1587 INTREPID2_TEST_FOR_EXCEPTION( inLeft.rank() != inRight.rank(), std::invalid_argument,
1588 ">>> ERROR (RealSpaceTools::vecprod): Right and left arrays must be have the same rank!");
1589 INTREPID2_TEST_FOR_EXCEPTION( inLeft.rank() != vecProd.rank(), std::invalid_argument,
1590 ">>> ERROR (RealSpaceTools::vecprod): Left and vecProd arrays must be have the same rank!");
1591 for (size_type i=0;i<inLeft.rank();++i) {
1592 INTREPID2_TEST_FOR_EXCEPTION( inLeft.extent(i) != inRight.extent(i), std::invalid_argument,
1593 ">>> ERROR (RealSpaceTools::vecprod): Dimensions of matrix arguments do not agree!");
1594 INTREPID2_TEST_FOR_EXCEPTION( inLeft.extent(i) != vecProd.extent(i), std::invalid_argument,
1595 ">>> ERROR (RealSpaceTools::vecprod): Dimensions of matrix arguments do not agree!");
1596 }
1597
1598 // (2) spatial dimension ordinal = array rank - 1. Suffices to check only one array because we just
1599 // checked whether or not the arrays have matching dimensions.
1600 INTREPID2_TEST_FOR_EXCEPTION( inLeft.extent(inLeft.rank()-1) < 2 ||
1601 inLeft.extent(inLeft.rank()-1) > 3, std::invalid_argument,
1602 ">>> ERROR (RealSpaceTools::vecprod): Dimensions of arrays (rank-1) must be 2 or 3!");
1603 }
1604#endif
1605 using MemSpaceType = typename DeviceType::memory_space;
1606 constexpr bool are_accessible =
1607 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1608 typename decltype(vecProd)::memory_space>::accessible &&
1609 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1610 typename decltype(inLeft)::memory_space>::accessible &&
1611 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1612 typename decltype(inRight)::memory_space>::accessible;
1613 static_assert(are_accessible, "RealSpaceTools<DeviceType>::vecprod(..): input/output views' memory spaces are not compatible with DeviceType");
1614
1615 using FunctorType = FunctorRealSpaceTools::F_vecprod<decltype(vecProd),decltype(inLeft),decltype(inRight)>;
1616
1617 const auto r = inLeft.rank();
1618 const auto loopSize = ( r == 1 ? 1 :
1619 r == 2 ? inLeft.extent(0) :
1620 inLeft.extent(0)*inLeft.extent(1) );
1621 const bool is_vecprod_3d = (inLeft.extent(inLeft.rank() - 1) == 3);
1622 using ExecSpaceType = typename DeviceType::execution_space;
1623 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
1624 Kokkos::parallel_for( policy, FunctorType(vecProd, inLeft, inRight, is_vecprod_3d) );
1625 }
1626
1627 // ------------------------------------------------------------------------------------
1628
1629} // namespace Intrepid2
1630
1631#endif
1632
1633
1634
1635
1636
1637// =====================================
1638// Too much things...
1639//
1640
1641// template<class ...inVec1Properties,
1642// class ...inVec2Properties>
1643// KOKKOS_INLINE_FUNCTION
1644// typename Kokkos::DynRankView<inVec1Properties...>::value_type
1645// RealSpaceTools<DeviceType>::
1646// dot( const Kokkos::DynRankView<inVec1Properties...> inVec1,
1647// const Kokkos::DynRankView<inVec2Properties...> inVec2 ) {
1648
1649// #ifdef HAVE_INTREPID2_DEBUG
1650// INTREPID2_TEST_FOR_ABORT( inVec1.rank != 1 || inVec2.rank() != 1,
1651// ">>> ERROR (RealSpaceTools::dot): Vector arguments must have rank 1!");
1652// INTREPID2_TEST_FOR_ABORT( inVec1.extent(0) != inVec2.extent(0),
1653// ">>> ERROR (RealSpaceTools::dot): Dimensions of vector arguments must agree!");
1654// #endif
1655// typedef typename Kokkos::DynRankView<inVec1Properties...>::value_type value_type;
1656
1657// // designed for small problems
1658// value_type r_val(0);
1659
1660// // * This is probably not necessary
1661// if (std::is_same<ExecSpace,Kokkos::Serial>::value) {
1662// const ordinal_type iend = iVec1.extent(0);
1663// for (ordinal_type i=0;i<iend;++i)
1664// r_val += inVec1(i)*inVec2(i);
1665// } else {
1666// struct Functor {
1667// Kokkos::DynRankView<inVec1Properties...> _inVec1;
1668// Kokkos::DynRankView<inVec2Properties...> _inVec2;
1669
1670// KOKKOS_INLINE_FUNCTION
1671// Functor(const Kokkos::DynRankView<inVec1Properties...> inVec1_,
1672// const Kokkos::DynRankView<inVec2Properties...> inVec2_);
1673
1674// KOKKOS_INLINE_FUNCTION
1675// void operator()(ordinal_type i, value_type &dst ) const {
1676// dst += _inVec1(i)*_inVec2(i);
1677// }
1678
1679// KOKKOS_INLINE_FUNCTION
1680// void join(volatile value_type &dst ,
1681// const volatile value_type &src) const {
1682// dst += src;
1683// }
1684// };
1685// const ordinal_type iend = inVec1.extent(0);
1686// Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, iend);
1687// Kokkos::parallel_for( policy, Functor(inVec1, inVec2), r_val );
1688// }
1689// return r_val;
1690// }
1691
1692
enable_if_t< has_rank_method< Functor >::value, unsigned > KOKKOS_INLINE_FUNCTION getFunctorRank(const Functor &functor)
KOKKOS_FORCEINLINE_FUNCTION constexpr std::enable_if<!std::is_pod< 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....
Implementation of basic linear algebra functionality in Euclidean space.
static void extractScalarValues(Kokkos::DynRankView< outputValueType, outputProperties... > output, const Kokkos::DynRankView< inputValueType, inputProperties... > input)
Extract scalar type values from Sacado-based array.
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.
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),...
static void inverse(InverseMatrixViewType inverseMats, MatrixViewType inMats)
Computes inverses of nonsingular matrices stored in an array of total rank 2 (single matrix),...
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.
static void absval(Kokkos::DynRankView< absArrayValueType, absArrayProperties... > absArray, const Kokkos::DynRankView< inArrayValueType, inArrayProperties... > inArray)
Computes absolute value of an array.
static void AtA(Kokkos::DynRankView< outMatValueType, outMatProperties... > outMats, const Kokkos::DynRankView< inMatValueType, inMatProperties... > inMats)
Computes the matrix-matrix product , for an input rectangular matrix .
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
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.
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.
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),...
static void det(DeterminantArrayViewType detArray, const MatrixViewType inMats)
Computes determinants of matrices stored in an array of total rank 3 (array of matrices),...
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),...
static void clone(Kokkos::DynRankView< outputValueType, outputProperties... > output, const Kokkos::DynRankView< inputValueType, inputProperties... > input)
Clone input array.
small utility functions
Functor to compute matvec see Intrepid2::RealSpaceTools for more.
Functor to compute absolute value see Intrepid2::RealSpaceTools for more.
Functor to add md arrays see Intrepid2::RealSpaceTools for more.
Functor for clone see Intrepid2::RealSpaceTools for more.
Functor to compute determinant see Intrepid2::RealSpaceTools for more.
Functor to compute dot product see Intrepid2::RealSpaceTools for more.
Functor for extractScalarValues see Intrepid2::RealSpaceTools for more.
Functor to compute inverse see Intrepid2::RealSpaceTools for more.
Functor to compute matvec see Intrepid2::RealSpaceTools for more.
Functor to scale md arrays see Intrepid2::RealSpaceTools for more.
Functor to subtract md arrays see Intrepid2::RealSpaceTools for more.
Functor to compute transpose see Intrepid2::RealSpaceTools for more.
Functor to compute vecprod see Intrepid2::RealSpaceTools for more.
Functor to compute vector norm see Intrepid2::RealSpaceTools for more.
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.
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.
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.