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