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 
52 namespace 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 
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 
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 
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>
468  struct F_vectorNorm {
469  normArrayViewType _normArray;
470  inVecViewType _inVecs;
471  const ENorm _normType;
472 
473  KOKKOS_INLINE_FUNCTION
474  F_vectorNorm( normArrayViewType normArray_,
475  inVecViewType inVecs_,
476  const ENorm normType_ )
477  : _normArray(normArray_), _inVecs(inVecs_), _normType(normType_) {}
478 
479  KOKKOS_INLINE_FUNCTION
480  void operator()(const ordinal_type iter) const {
481  ordinal_type i, j;
482  unrollIndex( i, j,
483  _normArray.extent(0),
484  _normArray.extent(1),
485  iter );
486 
487  auto vec = ( _inVecs.rank() == 2 ? Kokkos::subview(_inVecs, i, Kokkos::ALL()) :
488  Kokkos::subview(_inVecs, i, j, Kokkos::ALL()) );
489 
490  _normArray(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");
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");
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 
848 template<class DeterminantArrayViewType, class MatrixViewType>
849 static void
850 det( 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 
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 
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 
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 
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 
1332  template<typename DeviceType>
1333  template<typename matVecValueType, class ...matVecProperties,
1334  typename inMatValueType, class ...inMatProperties,
1335  typename inVecValueType, class ...inVecProperties>
1336  void
1338  matvec( Kokkos::DynRankView<matVecValueType,matVecProperties...> matVecs,
1339  const Kokkos::DynRankView<inMatValueType, inMatProperties...> inMats,
1340  const Kokkos::DynRankView<inVecValueType, inVecProperties...> inVecs ) {
1341 
1342 #ifdef HAVE_INTREPID2_DEBUG
1343  INTREPID2_TEST_FOR_EXCEPTION( inMats.rank() < 2 || inMats.rank() > 4, std::invalid_argument,
1344  ">>> ERROR (RealSpaceTools::matvec): Rank of matrix array must be 2, 3 or 4!");
1345  INTREPID2_TEST_FOR_EXCEPTION( matVecs.rank() < 1 || matVecs.rank() > 3, std::invalid_argument,
1346  ">>> ERROR (RealSpaceTools::matvec): Rank of matVecs array must be 1, 2 or 3!");
1347  INTREPID2_TEST_FOR_EXCEPTION( inVecs.rank() < 1 || inVecs.rank() > 3, std::invalid_argument,
1348  ">>> ERROR (RealSpaceTools::matvec): Rank of inVecs array must be 1, 2 or 3!");
1349  if (inMats.rank() == 2) {
1350  // a single matrix, multiple input and output
1351  INTREPID2_TEST_FOR_EXCEPTION( matVecs.rank() != inVecs.rank(), std::invalid_argument,
1352  ">>> ERROR (RealSpaceTools::matvec): Vector arrays must be have the same rank!");
1353  // output must match
1354  for (ordinal_type i=0;i< (static_cast<ordinal_type>(inVecs.rank())-1);++i) {
1355  INTREPID2_TEST_FOR_EXCEPTION( matVecs.extent(i) != inVecs.extent(i), std::invalid_argument,
1356  ">>> ERROR (RealSpaceTools::matvec): Dimensions of vector array arguments do not agree!");
1357  }
1358  // matvec compatibility
1359  INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(0) != matVecs.extent(matVecs.rank()-1) ||
1360  inMats.extent(1) != inVecs.extent(inVecs.rank()-1), std::invalid_argument,
1361  ">>> ERROR (RealSpaceTools::matvec): Matvec dimensions are not compatible each other.");
1362  } else if (inVecs.rank() == 1) {
1363  // multiple matrix, single input and multiple output
1364  INTREPID2_TEST_FOR_EXCEPTION( inMats.rank() != (matVecs.rank()+1), std::invalid_argument,
1365  ">>> ERROR (RealSpaceTools::matvec): The result vector and matrix array arguments do not have compatible ranks!");
1366  for (ordinal_type i=0;i<(static_cast<ordinal_type>(inMats.rank())-2);++i) {
1367  INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(i) != matVecs.extent(i), std::invalid_argument,
1368  ">>> ERROR (RealSpaceTools::matvec): Dimensions of vector and matrix array arguments do not agree!");
1369  }
1370  // matvec compatibility
1371  INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(inMats.rank()-2) != matVecs.extent(matVecs.rank()-1) ||
1372  inMats.extent(inMats.rank()-1) != inVecs.extent(0), std::invalid_argument,
1373  ">>> ERROR (RealSpaceTools::matvec): Matvec dimensions are not compatible each other.");
1374  } else {
1375  // multiple matrix, multiple input and multiple output
1376  INTREPID2_TEST_FOR_EXCEPTION( inMats.rank() != (matVecs.rank()+1), std::invalid_argument,
1377  ">>> ERROR (RealSpaceTools::matvec): The result vector and matrix array arguments do not have compatible ranks!");
1378  INTREPID2_TEST_FOR_EXCEPTION( matVecs.rank() != inVecs.rank(), std::invalid_argument,
1379  ">>> ERROR (RealSpaceTools::matvec): Vector arrays must be have the same rank!");
1380  for (ordinal_type i=0;i<(static_cast<ordinal_type>(inMats.rank())-2);++i) {
1381  INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(i) != matVecs.extent(i), std::invalid_argument,
1382  ">>> ERROR (RealSpaceTools::matvec): Dimensions of vector and matrix array arguments do not agree!");
1383  }
1384  for (size_type i=0;i<inVecs.rank();++i) {
1385  INTREPID2_TEST_FOR_EXCEPTION( matVecs.extent(i) != inVecs.extent(i), std::invalid_argument,
1386  ">>> ERROR (RealSpaceTools::matvec): Dimensions of vector array arguments do not agree!");
1387  }
1388  INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(inMats.rank()-1) != inVecs.extent(inVecs.rank()-1), std::invalid_argument,
1389  ">>> ERROR (RealSpaceTools::matvec): Matrix column dimension does not match to the length of a vector!");
1390  }
1391 #endif
1392  using MemSpaceType = typename DeviceType::memory_space;
1393  constexpr bool are_accessible =
1394  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1395  typename decltype(matVecs)::memory_space>::accessible &&
1396  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1397  typename decltype(inMats)::memory_space>::accessible &&
1398  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1399  typename decltype(inVecs)::memory_space>::accessible;
1400  static_assert(are_accessible, "RealSpaceTools<DeviceType>::matvec(..): input/output views' memory spaces are not compatible with DeviceType");
1401 
1403 
1404  size_type loopSize = 1;
1405  const ordinal_type r = matVecs.rank() - 1;
1406  for (ordinal_type i=0;i<r;++i)
1407  loopSize *= matVecs.extent(i);
1408 
1409  using ExecSpaceType = typename DeviceType::execution_space;
1410  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
1411  Kokkos::parallel_for( policy, FunctorType(matVecs, inMats, inVecs) );
1412  }
1413 
1414  // ------------------------------------------------------------------------------------
1415 
1416  namespace FunctorRealSpaceTools {
1420  template<typename vecProdViewType,
1421  typename inLeftViewType,
1422  typename inRightViewType>
1423  struct F_vecprod {
1424  vecProdViewType _vecProd;
1425  const inLeftViewType _inLeft;
1426  const inRightViewType _inRight;
1427  const bool _is_vecprod_3d;
1428 
1429  KOKKOS_INLINE_FUNCTION
1430  F_vecprod( vecProdViewType vecProd_,
1431  inLeftViewType inLeft_,
1432  inRightViewType inRight_,
1433  const bool is_vecprod_3d_ )
1434  : _vecProd(vecProd_), _inLeft(inLeft_), _inRight(inRight_), _is_vecprod_3d(is_vecprod_3d_) {}
1435 
1436  KOKKOS_INLINE_FUNCTION
1437  void operator()(const size_type iter) const {
1438  ordinal_type i, j;
1439  unrollIndex( i, j,
1440  _inLeft.extent(0),
1441  _inLeft.extent(1),
1442  iter );
1443 
1444  const auto r = _inLeft.rank();
1445 
1446  auto left = ( r == 1 ? Kokkos::subview(_inLeft, Kokkos::ALL()) :
1447  r == 2 ? Kokkos::subview(_inLeft, i, Kokkos::ALL()) :
1448  Kokkos::subview(_inLeft, i, j, Kokkos::ALL()) );
1449 
1450  auto right = ( r == 1 ? Kokkos::subview(_inRight, Kokkos::ALL()) :
1451  r == 2 ? Kokkos::subview(_inRight, i, Kokkos::ALL()) :
1452  Kokkos::subview(_inRight, i, j, Kokkos::ALL()) );
1453 
1454  auto result = ( r == 1 ? Kokkos::subview(_vecProd, Kokkos::ALL()) :
1455  r == 2 ? Kokkos::subview(_vecProd, i, Kokkos::ALL()) :
1456  Kokkos::subview(_vecProd, i, j, Kokkos::ALL()) );
1457 
1458  // branch prediction is possible
1459  if (_is_vecprod_3d) {
1460  result(0) = left(1)*right(2) - left(2)*right(1);
1461  result(1) = left(2)*right(0) - left(0)*right(2);
1462  result(2) = left(0)*right(1) - left(1)*right(0);
1463  } else {
1464  result(0) = left(0)*right(1) - left(1)*right(0);
1465  }
1466  }
1467  };
1468  }
1469 
1470  template<typename DeviceType>
1471  template<typename vecProdValueType, class ...vecProdProperties,
1472  typename inLeftValueType, class ...inLeftProperties,
1473  typename inRightValueType, class ...inRightProperties>
1474  void
1476  vecprod( Kokkos::DynRankView<vecProdValueType,vecProdProperties...> vecProd,
1477  const Kokkos::DynRankView<inLeftValueType, inLeftProperties...> inLeft,
1478  const Kokkos::DynRankView<inRightValueType,inRightProperties...> inRight ) {
1479 
1480 #ifdef HAVE_INTREPID2_DEBUG
1481  {
1482  /*
1483  * Check array rank and spatial dimension range (if applicable)
1484  * (1) all array arguments are required to have matching dimensions and rank: (D), (I0,D) or (I0,I1,D)
1485  * (2) spatial dimension should be 2 or 3
1486  */
1487  // (1) check rank range on inLeft and then compare the other arrays with inLeft
1488  INTREPID2_TEST_FOR_EXCEPTION( inLeft.rank() < 1 || inLeft.rank() > 3, std::invalid_argument,
1489  ">>> ERROR (RealSpaceTools::vecprod): Rank of matrix array must be 1, 2, or 3!");
1490  INTREPID2_TEST_FOR_EXCEPTION( inLeft.rank() != inRight.rank(), std::invalid_argument,
1491  ">>> ERROR (RealSpaceTools::vecprod): Right and left arrays must be have the same rank!");
1492  INTREPID2_TEST_FOR_EXCEPTION( inLeft.rank() != vecProd.rank(), std::invalid_argument,
1493  ">>> ERROR (RealSpaceTools::vecprod): Left and vecProd arrays must be have the same rank!");
1494  for (size_type i=0;i<inLeft.rank();++i) {
1495  INTREPID2_TEST_FOR_EXCEPTION( inLeft.extent(i) != inRight.extent(i), std::invalid_argument,
1496  ">>> ERROR (RealSpaceTools::vecprod): Dimensions of matrix arguments do not agree!");
1497  INTREPID2_TEST_FOR_EXCEPTION( inLeft.extent(i) != vecProd.extent(i), std::invalid_argument,
1498  ">>> ERROR (RealSpaceTools::vecprod): Dimensions of matrix arguments do not agree!");
1499  }
1500 
1501  // (2) spatial dimension ordinal = array rank - 1. Suffices to check only one array because we just
1502  // checked whether or not the arrays have matching dimensions.
1503  INTREPID2_TEST_FOR_EXCEPTION( inLeft.extent(inLeft.rank()-1) < 2 ||
1504  inLeft.extent(inLeft.rank()-1) > 3, std::invalid_argument,
1505  ">>> ERROR (RealSpaceTools::vecprod): Dimensions of arrays (rank-1) must be 2 or 3!");
1506  }
1507 #endif
1508  using MemSpaceType = typename DeviceType::memory_space;
1509  constexpr bool are_accessible =
1510  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1511  typename decltype(vecProd)::memory_space>::accessible &&
1512  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1513  typename decltype(inLeft)::memory_space>::accessible &&
1514  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1515  typename decltype(inRight)::memory_space>::accessible;
1516  static_assert(are_accessible, "RealSpaceTools<DeviceType>::vecprod(..): input/output views' memory spaces are not compatible with DeviceType");
1517 
1519 
1520  const auto r = inLeft.rank();
1521  const auto loopSize = ( r == 1 ? 1 :
1522  r == 2 ? inLeft.extent(0) :
1523  inLeft.extent(0)*inLeft.extent(1) );
1524  const bool is_vecprod_3d = (inLeft.extent(inLeft.rank() - 1) == 3);
1525  using ExecSpaceType = typename DeviceType::execution_space;
1526  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
1527  Kokkos::parallel_for( policy, FunctorType(vecProd, inLeft, inRight, is_vecprod_3d) );
1528  }
1529 
1530  // ------------------------------------------------------------------------------------
1531 
1532 } // namespace Intrepid2
1533 
1534 #endif
1535 
1536 
1537 
1538 
1539 
1540 // =====================================
1541 // Too much things...
1542 //
1543 
1544 // template<class ...inVec1Properties,
1545 // class ...inVec2Properties>
1546 // KOKKOS_INLINE_FUNCTION
1547 // typename Kokkos::DynRankView<inVec1Properties...>::value_type
1548 // RealSpaceTools<DeviceType>::
1549 // dot( const Kokkos::DynRankView<inVec1Properties...> inVec1,
1550 // const Kokkos::DynRankView<inVec2Properties...> inVec2 ) {
1551 
1552 // #ifdef HAVE_INTREPID2_DEBUG
1553 // INTREPID2_TEST_FOR_ABORT( inVec1.rank != 1 || inVec2.rank() != 1,
1554 // ">>> ERROR (RealSpaceTools::dot): Vector arguments must have rank 1!");
1555 // INTREPID2_TEST_FOR_ABORT( inVec1.extent(0) != inVec2.extent(0),
1556 // ">>> ERROR (RealSpaceTools::dot): Dimensions of vector arguments must agree!");
1557 // #endif
1558 // typedef typename Kokkos::DynRankView<inVec1Properties...>::value_type value_type;
1559 
1560 // // designed for small problems
1561 // value_type r_val(0);
1562 
1563 // // * This is probably not necessary
1564 // if (std::is_same<ExecSpace,Kokkos::Serial>::value) {
1565 // const ordinal_type iend = iVec1.extent(0);
1566 // for (ordinal_type i=0;i<iend;++i)
1567 // r_val += inVec1(i)*inVec2(i);
1568 // } else {
1569 // struct Functor {
1570 // Kokkos::DynRankView<inVec1Properties...> _inVec1;
1571 // Kokkos::DynRankView<inVec2Properties...> _inVec2;
1572 
1573 // KOKKOS_INLINE_FUNCTION
1574 // Functor(const Kokkos::DynRankView<inVec1Properties...> inVec1_,
1575 // const Kokkos::DynRankView<inVec2Properties...> inVec2_);
1576 
1577 // KOKKOS_INLINE_FUNCTION
1578 // void operator()(ordinal_type i, value_type &dst ) const {
1579 // dst += _inVec1(i)*_inVec2(i);
1580 // }
1581 
1582 // KOKKOS_INLINE_FUNCTION
1583 // void join(volatile value_type &dst ,
1584 // const volatile value_type &src) const {
1585 // dst += src;
1586 // }
1587 // };
1588 // const ordinal_type iend = inVec1.extent(0);
1589 // Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, iend);
1590 // Kokkos::parallel_for( policy, Functor(inVec1, inVec2), r_val );
1591 // }
1592 // return r_val;
1593 // }
1594 
1595 
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 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).
small utility functions
Functor to compute dot product see Intrepid2::RealSpaceTools for more.
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 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 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 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.
Functor to compute inverse see Intrepid2::RealSpaceTools for more.
Functor to compute determinant see Intrepid2::RealSpaceTools for more.
ENorm
Enumeration of norm types for vectors and functions.
static void inverse(InverseMatrixViewType inverseMats, MatrixViewType inMats)
Computes inverses of nonsingular matrices stored in an array of total rank 2 (single matrix)...
Functor to compute vecprod see Intrepid2::RealSpaceTools for more.
Functor to compute absolute value see Intrepid2::RealSpaceTools for more.
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).
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.
Functor to subtract md arrays 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.
Functor to scale md arrays see Intrepid2::RealSpaceTools for more.
Functor to add md arrays see Intrepid2::RealSpaceTools for more.
Functor to compute vector norm see Intrepid2::RealSpaceTools for more.
Functor for extractScalarValues see Intrepid2::RealSpaceTools for more.
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 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
Functor to compute transpose see Intrepid2::RealSpaceTools for more.
KOKKOS_FORCEINLINE_FUNCTION constexpr std::enable_if< !std::is_pod< T >::value, typename ScalarTraits< T >::scalar_type >::type get_scalar_value(const T &obj)
functions returning the scalar value. for pod types, they return the input object itself...
static void extractScalarValues(Kokkos::DynRankView< outputValueType, outputProperties... > output, const Kokkos::DynRankView< inputValueType, inputProperties... > input)
Extract scalar type values from Sacado-based array.
enable_if_t< has_rank_method< Functor >::value, unsigned > KOKKOS_INLINE_FUNCTION getFunctorRank(const Functor &functor)
Functor for clone see Intrepid2::RealSpaceTools for more.
Functor to compute matvec see Intrepid2::RealSpaceTools for more.
static void clone(Kokkos::DynRankView< outputValueType, outputProperties... > output, const Kokkos::DynRankView< inputValueType, inputProperties... > input)
Clone input 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.