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
66  RealSpaceTools<DeviceType>::Serial::
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
136  RealSpaceTools<DeviceType>::Serial::
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
184  RealSpaceTools<DeviceType>::Serial::
185  dot( const Kokkos::DynRankView<inVec1ValueType,inVec1Properties...> inVec1,
186  const Kokkos::DynRankView<inVec2ValueType,inVec2Properties...> inVec2 ) {
187 
188 #ifdef HAVE_INTREPID2_DEBUG
189  {
190  bool dbgInfo = false;
191  INTREPID2_TEST_FOR_DEBUG_ABORT( inVec1.rank() != 1 || inVec2.rank() != 1, dbgInfo,
192  ">>> ERROR (RealSpaceTools::dot): Vector arguments must have rank 1!");
193  INTREPID2_TEST_FOR_DEBUG_ABORT( inVec1.extent(0) != inVec2.extent(0), dbgInfo,
194  ">>> ERROR (RealSpaceTools::dot): Dimensions of vector arguments must agree!");
195 #ifdef INTREPID2_TEST_FOR_EXCEPTION_OVERRIDE_TO_CONTINUE
196  if (dbgInfo) return inVec1ValueType(0);
197 #endif
198  }
199 #endif
200  typedef inVec1ValueType value_type;
201 
202  // designed for small problems
203  value_type r_val(0);
204 
205  const ordinal_type iend = inVec1.extent(0);
206  for (ordinal_type i=0;i<iend;++i)
207  r_val += inVec1(i)*inVec2(i);
208 
209  return r_val;
210  }
211 
212  // ------------------------------------------------------------------------------------
213 
214  //
215  // use parallel for
216  //
217 
218  // ------------------------------------------------------------------------------------
219 
220  namespace FunctorRealSpaceTools {
221 
225  template<typename OutputViewType,
226  typename inputViewType>
227  struct F_extractScalarValues {
228  OutputViewType _output;
229  inputViewType _input;
230 
231  KOKKOS_INLINE_FUNCTION
232  F_extractScalarValues( OutputViewType output_,
233  inputViewType input_ )
234  : _output(output_), _input(input_) {}
235 
236  KOKKOS_INLINE_FUNCTION
237  void operator()(const ordinal_type i) const {
238  const ordinal_type jend = _input.extent(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
256  RealSpaceTools<DeviceType>::
257  extractScalarValues( Kokkos::DynRankView<outputValueType,outputProperties...> output,
258  const Kokkos::DynRankView<inputValueType, inputProperties...> input ) {
259 
260  using MemSpaceType = typename DeviceType::memory_space;
261  constexpr bool are_accessible =
262  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
263  typename decltype(output)::memory_space>::accessible &&
264  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
265  typename decltype(input)::memory_space>::accessible;
266  static_assert(are_accessible, "RealSpaceTools<DeviceType>::extractScalarValues(..): input/output views' memory spaces are not compatible with DeviceType");
267 
268  using FunctorType = FunctorRealSpaceTools::F_extractScalarValues<decltype(output),decltype(input)>;
269 
270  const auto loopSize = input.extent(0);
271  using ExecSpaceType = typename DeviceType::execution_space;
272  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
273  Kokkos::parallel_for( policy, FunctorType(output, input) );
274  }
275 
276  namespace FunctorRealSpaceTools {
280  template<typename OutputViewType,
281  typename inputViewType>
282  struct F_clone {
283  OutputViewType _output;
284  inputViewType _input;
285 
286  KOKKOS_INLINE_FUNCTION
287  F_clone( OutputViewType output_,
288  inputViewType input_ )
289  : _output(output_), _input(input_) {}
290 
291  KOKKOS_INLINE_FUNCTION
292  void operator()(const ordinal_type iter) const {
293  ordinal_type k0(0), k1(0), k2(0);
294  const auto rankDiff = _output.rank() - _input.rank();
295  switch (rankDiff) {
296  case 3:
297  unrollIndex( k0, k1, k2,
298  _output.extent(0),
299  _output.extent(1),
300  _output.extent(2),
301  iter );
302  break;
303  case 2:
304  unrollIndex( k0, k1,
305  _output.extent(0),
306  _output.extent(1),
307  iter );
308  break;
309  case 1:
310  k0 = iter;
311  break;
312  }
313 
314  auto out = (rankDiff == 3 ? Kokkos::subview(_output, k0, k1, k2, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()) :
315  rankDiff == 2 ? Kokkos::subview(_output, k0, k1, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()) :
316  rankDiff == 1 ? Kokkos::subview(_output, k0, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()) :
317  Kokkos::subview(_output, Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()));
318  const ordinal_type iend = _input.extent(0);
319  const ordinal_type jend = _input.extent(1);
320  const ordinal_type kend = _input.extent(2);
321  for (ordinal_type i=0;i<iend;++i)
322  for (ordinal_type j=0;j<jend;++j)
323  for (ordinal_type k=0;k<kend;++k)
324  out.access(i,j,k) = _input.access(i,j,k);
325  }
326  };
327  }
328 
329  template<typename DeviceType>
330  template<typename outputValueType, class ...outputProperties,
331  typename inputValueType, class ...inputProperties>
332  void
333  RealSpaceTools<DeviceType>::
334  clone( Kokkos::DynRankView<outputValueType,outputProperties...> output,
335  const Kokkos::DynRankView<inputValueType,inputProperties...> input ) {
336 #ifdef HAVE_INTREPID2_DEBUG
337  {
338  // input has at most rank 3
339  INTREPID2_TEST_FOR_EXCEPTION( input.rank() > 3, std::invalid_argument,
340  ">>> ERROR (RealSpaceTools::clone): Input container has rank larger than 3.");
341 
342 
343  INTREPID2_TEST_FOR_EXCEPTION( input.rank() > output.rank(), std::invalid_argument,
344  ">>> ERROR (RealSpaceTools::clone): Input container rank should be smaller than ouput rank.");
345 
346  const size_type rankDiff = output.rank() - input.rank();
347 
348  // Difference of output and input rank is at most 3.
349  INTREPID2_TEST_FOR_EXCEPTION( rankDiff > 3, std::invalid_argument,
350  ">>> ERROR (RealSpaceTools::clone): Difference between output container rank and input container rank is larger than 3.");
351 
352 
353  for (size_type i=0;i<input.rank();++i) {
354  INTREPID2_TEST_FOR_EXCEPTION( input.extent(i) != output.extent(rankDiff + i), std::invalid_argument,
355  ">>> ERROR (RealSpaceTools::clone): Dimensions of array arguments do not agree!");
356  }
357  }
358 #endif
359  using MemSpaceType = typename DeviceType::memory_space;
360  constexpr bool are_accessible =
361  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
362  typename decltype(output)::memory_space>::accessible &&
363  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
364  typename decltype(input)::memory_space>::accessible;
365  static_assert(are_accessible, "RealSpaceTools<DeviceType>::clone(..): input/output views' memory spaces are not compatible with DeviceType");
366 
367  using FunctorType = FunctorRealSpaceTools::F_clone<decltype(output),decltype(input)>;
368 
369  size_type loopSize = 1;
370  const ordinal_type out_rank = output.rank();
371  const ordinal_type in_rank = input.rank();
372  const ordinal_type rankDiff = out_rank - in_rank;
373  for (ordinal_type i=0;i<rankDiff;++i)
374  loopSize *= output.extent(i);
375 
376  using ExecSpaceType = typename DeviceType::execution_space;
377  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
378  Kokkos::parallel_for( policy, FunctorType(output, input) );
379  }
380 
381  namespace FunctorRealSpaceTools {
385  template<typename absArrayViewType,
386  typename inArrayViewType>
387  struct F_absval {
388  absArrayViewType _absArray;
389  inArrayViewType _inArray;
390 
391  KOKKOS_INLINE_FUNCTION
392  F_absval( absArrayViewType absArray_,
393  inArrayViewType inArray_ )
394  : _absArray(absArray_), _inArray(inArray_) {}
395 
396  KOKKOS_INLINE_FUNCTION
397  void operator()(const ordinal_type i) const {
398  const ordinal_type jend = _inArray.extent(1);
399  const ordinal_type kend = _inArray.extent(2);
400  const ordinal_type lend = _inArray.extent(3);
401  const ordinal_type mend = _inArray.extent(4);
402 
403  typedef typename inArrayViewType::non_const_value_type value_type;
404 
405  for (ordinal_type j=0;j<jend;++j)
406  for (ordinal_type k=0;k<kend;++k)
407  for (ordinal_type l=0;l<lend;++l)
408  for (ordinal_type m=0;m<mend;++m)
409  _absArray.access(i,j,k,l,m) = Util<value_type>::abs(_inArray.access(i,j,k,l,m));
410  }
411  };
412  }
413 
414  template<typename DeviceType>
415  template<typename absArrayValueType, class ...absArrayProperties,
416  typename inArrayValueType, class ...inArrayProperties>
417  void
418  RealSpaceTools<DeviceType>::
419  absval( Kokkos::DynRankView<absArrayValueType,absArrayProperties...> absArray,
420  const Kokkos::DynRankView<inArrayValueType, inArrayProperties...> inArray ) {
421 #ifdef HAVE_INTREPID2_DEBUG
422  {
423  INTREPID2_TEST_FOR_EXCEPTION( inArray.rank() > 5, std::invalid_argument,
424  ">>> ERROR (RealSpaceTools::absval): Input array container has rank larger than 5.");
425 
426  INTREPID2_TEST_FOR_EXCEPTION( inArray.rank() != absArray.rank(), std::invalid_argument,
427  ">>> ERROR (RealSpaceTools::absval): Array arguments must have identical ranks!");
428  for (size_type i=0;i<inArray.rank();++i) {
429  INTREPID2_TEST_FOR_EXCEPTION( inArray.extent(i) != absArray.extent(i), std::invalid_argument,
430  ">>> ERROR (RealSpaceTools::absval): Dimensions of array arguments do not agree!");
431  }
432  }
433 #endif
434  using MemSpaceType = typename DeviceType::memory_space;
435  constexpr bool are_accessible =
436  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
437  typename decltype(absArray)::memory_space>::accessible &&
438  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
439  typename decltype(inArray)::memory_space>::accessible;
440  static_assert(are_accessible, "RealSpaceTools<DeviceType>::absval(..): input/output views' memory spaces are not compatible with DeviceType");
441 
442  using FunctorType = FunctorRealSpaceTools::F_absval<decltype(absArray),decltype(inArray)>;
443 
444  const auto loopSize = inArray.extent(0);
445  using ExecSpaceType = typename DeviceType::execution_space;
446  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
447  Kokkos::parallel_for( policy, FunctorType(absArray, inArray) );
448  }
449 
450  // ------------------------------------------------------------------------------------
451 
452  template<typename DeviceType>
453  template<typename inoutArrayValueType, class ...inoutAbsArrayProperties>
454  void
455  RealSpaceTools<DeviceType>::
456  absval( Kokkos::DynRankView<inoutArrayValueType,inoutAbsArrayProperties...> inoutAbsArray ) {
457  RealSpaceTools<DeviceType>::absval(inoutAbsArray, inoutAbsArray);
458  }
459 
460  // ------------------------------------------------------------------------------------
461 
462  namespace FunctorRealSpaceTools {
466  template<typename normArrayViewType,
467  typename inVecViewType>
468  struct F_vectorNorm {
469  normArrayViewType _normArray;
470  inVecViewType _inVecs;
471  const ENorm _normType;
472 
473  KOKKOS_INLINE_FUNCTION
474  F_vectorNorm( normArrayViewType normArray_,
475  inVecViewType inVecs_,
476  const ENorm normType_ )
477  : _normArray(normArray_), _inVecs(inVecs_), _normType(normType_) {}
478 
479  KOKKOS_INLINE_FUNCTION
480  void operator()(const ordinal_type iter) const {
481  ordinal_type i, j;
482  unrollIndex( i, j,
483  _normArray.extent(0),
484  _normArray.extent(1),
485  iter );
486 
487  auto vec = ( _inVecs.rank() == 2 ? Kokkos::subview(_inVecs, i, Kokkos::ALL()) :
488  Kokkos::subview(_inVecs, i, j, Kokkos::ALL()) );
489 
490  _normArray(i, j) = RealSpaceTools<>::Serial::vectorNorm(vec, _normType);
491  }
492  };
493  }
494 
495  template<typename DeviceType>
496  template<typename normArrayValueType, class ...normArrayProperties,
497  typename inVecValueType, class ...inVecProperties>
498  void
499  RealSpaceTools<DeviceType>::
500  vectorNorm( Kokkos::DynRankView<normArrayValueType,normArrayProperties...> normArray,
501  const Kokkos::DynRankView<inVecValueType, inVecProperties...> inVecs,
502  const ENorm normType ) {
503 #ifdef HAVE_INTREPID2_DEBUG
504  {
505  INTREPID2_TEST_FOR_EXCEPTION( inVecs.rank() != (normArray.rank()+1), std::invalid_argument,
506  ">>> ERROR (RealSpaceTools::vectorNorm): Ranks of norm and vector array arguments are incompatible!");
507  INTREPID2_TEST_FOR_EXCEPTION( inVecs.rank() < 2 || inVecs.rank() > 3, std::invalid_argument,
508  ">>> ERROR (RealSpaceTools::vectorNorm): Rank of vector array must be 2 or 3!");
509  for (size_type i=0;i<inVecs.rank()-1;++i) {
510  INTREPID2_TEST_FOR_EXCEPTION( inVecs.extent(i) != normArray.extent(i), std::invalid_argument,
511  ">>> ERROR (RealSpaceTools::vectorNorm): Dimensions of norm and vector arguments do not agree!");
512  }
513  }
514 #endif
515 
516  using MemSpaceType = typename DeviceType::memory_space;
517  constexpr bool are_accessible =
518  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
519  typename decltype(normArray)::memory_space>::accessible &&
520  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
521  typename decltype(inVecs)::memory_space>::accessible;
522  static_assert(are_accessible, "RealSpaceTools<DeviceType>::vectorNorm(..): input/output views' memory spaces are not compatible with DeviceType");
523  using FunctorType = FunctorRealSpaceTools::F_vectorNorm<decltype(normArray),decltype(inVecs)>;
524 
525  // normArray rank is either 1 or 2
526  const auto loopSize = normArray.extent(0)*normArray.extent(1);
527  using ExecSpaceType = typename DeviceType::execution_space;
528  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
529  Kokkos::parallel_for( policy, FunctorType(normArray, inVecs, normType) );
530  }
531 
532  // ------------------------------------------------------------------------------------
533 
534  namespace FunctorRealSpaceTools {
538  template<typename transposeMatViewType,
539  typename inMatViewType>
540  struct F_transpose {
541  transposeMatViewType _transposeMats;
542  inMatViewType _inMats;
543 
544  KOKKOS_INLINE_FUNCTION
545  F_transpose( transposeMatViewType transposeMats_,
546  inMatViewType inMats_)
547  : _transposeMats(transposeMats_), _inMats(inMats_) {}
548 
549  KOKKOS_INLINE_FUNCTION
550  void operator()(const size_type iter) const {
551  // the rank of normArray is either 1 or 2
552  const auto r = _transposeMats.rank();
553  ordinal_type _i = iter, _j = 0;
554 
555  if ( r > 3 )
556  unrollIndex( _i, _j,
557  _transposeMats.extent(0),
558  _transposeMats.extent(1),
559  iter );
560 
561  auto dst = ( r == 2 ? Kokkos::subview(_transposeMats, Kokkos::ALL(), Kokkos::ALL()) :
562  r == 3 ? Kokkos::subview(_transposeMats, _i, Kokkos::ALL(), Kokkos::ALL()) :
563  Kokkos::subview(_transposeMats, _i, _j, Kokkos::ALL(), Kokkos::ALL()) );
564 
565  auto src = ( r == 2 ? Kokkos::subview(_inMats, Kokkos::ALL(), Kokkos::ALL()) :
566  r == 3 ? Kokkos::subview(_inMats, _i, Kokkos::ALL(), Kokkos::ALL()) :
567  Kokkos::subview(_inMats, _i, _j, Kokkos::ALL(), Kokkos::ALL()) );
568 
569  for (size_type i=0;i<src.extent(0);++i) {
570  dst(i, i) = src(i, i);
571  for (size_type j=i+1;j<src.extent(1);++j) {
572  dst(i, j) = src(j, i);
573  dst(j, i) = src(i, j);
574  }
575  }
576  }
577  };
578  }
579 
580  template<typename DeviceType>
581  template<typename transposeMatValueType, class ...transposeMatProperties,
582  typename inMatValueType, class ...inMatProperties>
583  void
584  RealSpaceTools<DeviceType>::
585  transpose( Kokkos::DynRankView<transposeMatValueType,transposeMatProperties...> transposeMats,
586  const Kokkos::DynRankView<inMatValueType, inMatProperties...> inMats ) {
587 
588 #ifdef HAVE_INTREPID2_DEBUG
589  {
590 
591  INTREPID2_TEST_FOR_EXCEPTION( inMats.rank() != transposeMats.rank(), std::invalid_argument,
592  ">>> ERROR (RealSpaceTools::transpose): Matrix array arguments do not have identical ranks!");
593  INTREPID2_TEST_FOR_EXCEPTION( inMats.rank() < 2 || inMats.rank() > 4, std::invalid_argument,
594  ">>> ERROR (RealSpaceTools::transpose): Rank of matrix array must be 2, 3, or 4!");
595  for (size_type i=0;i<inMats.rank();++i) {
596  INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(i) != transposeMats.extent(i), std::invalid_argument,
597  ">>> ERROR (RealSpaceTools::transpose): Dimensions of matrix arguments do not agree!");
598  }
599  INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(inMats.rank()-2) != inMats.extent(inMats.rank()-1), std::invalid_argument,
600  ">>> ERROR (RealSpaceTools::transpose): Matrices are not square!");
601  }
602 #endif
603 
604  using MemSpaceType = typename DeviceType::memory_space;
605  constexpr bool are_accessible =
606  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
607  typename decltype(transposeMats)::memory_space>::accessible &&
608  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
609  typename decltype(inMats)::memory_space>::accessible;
610  static_assert(are_accessible, "RealSpaceTools<DeviceType>::transpose(..): input/output views' memory spaces are not compatible with DeviceType");
611  using FunctorType = FunctorRealSpaceTools::F_transpose<decltype(transposeMats),decltype(inMats)>;
612 
613  const auto r = transposeMats.rank();
614  const auto loopSize = ( r == 2 ? 1 :
615  r == 3 ? transposeMats.extent(0) :
616  transposeMats.extent(0)*transposeMats.extent(1) );
617 
618  using ExecSpaceType = typename DeviceType::execution_space;
619  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
620  Kokkos::parallel_for( policy, FunctorType(transposeMats, inMats) );
621  }
622 
623  // ------------------------------------------------------------------------------------
624 
625  namespace FunctorRealSpaceTools {
629  template<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
749  RealSpaceTools<DeviceType>::
750  inverse( InverseMatrixViewType inverseMats, MatrixViewType inMats ) {
751  const unsigned rank = getFunctorRank(inMats);
752 #ifdef HAVE_INTREPID2_DEBUG
753  {
754  INTREPID2_TEST_FOR_EXCEPTION( rank != getFunctorRank(inverseMats), std::invalid_argument,
755  ">>> ERROR (RealSpaceTools::inverse): Matrix array arguments do not have identical ranks!");
756  INTREPID2_TEST_FOR_EXCEPTION( rank < 3 || rank > 4, std::invalid_argument,
757  ">>> ERROR (RealSpaceTools::inverse): Rank of matrix array must be 3, or 4!");
758  for (size_type i=0;i<rank;++i) {
759  INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(i) != inverseMats.extent(i), std::invalid_argument,
760  ">>> ERROR (RealSpaceTools::inverse): Dimensions of matrix arguments do not agree!");
761  }
762  INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(rank-2) != inMats.extent(rank-1), std::invalid_argument,
763  ">>> ERROR (RealSpaceTools::inverse): Matrices are not square!");
764  INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(rank-2) < 1 || inMats.extent(rank-2) > 3, std::invalid_argument,
765  ">>> ERROR (RealSpaceTools::inverse): Spatial dimension must be 1, 2, or 3!");
766  }
767 #endif
768 
769  using ExecSpaceType = typename DeviceType::execution_space;
770  using FunctorType = FunctorRealSpaceTools::F_inverse<InverseMatrixViewType, MatrixViewType>;
771 
772  switch (rank) {
773  case 3: { // output P,D,D and input P,D,D
774  using range_policy_type = Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> >;
775  range_policy_type policy(0, inverseMats.extent(0));
776  Kokkos::parallel_for( policy, FunctorType(inverseMats, inMats) );
777  break;
778  }
779  case 4: { // output C,P,D,D and input C,P,D,D
780  using range_policy_type = Kokkos::MDRangePolicy
781  < ExecSpaceType, Kokkos::Rank<2>, Kokkos::IndexType<ordinal_type> >;
782  range_policy_type policy( { 0, 0 },
783  { inverseMats.extent(0), inverseMats.extent(1) } );
784  Kokkos::parallel_for( policy, FunctorType(inverseMats, inMats) );
785  break;
786  }
787  default: {
788  INTREPID2_TEST_FOR_EXCEPTION( true, std::invalid_argument,
789  ">>> ERROR (RealSpaceTools::inverse): Rank of matrix array must be 2, 3, or 4!");
790  }
791  }
792  }
793 
794  // ------------------------------------------------------------------------------------
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
855  RealSpaceTools<DeviceType>::
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
883  using FunctorType = FunctorRealSpaceTools::F_det<DeterminantArrayViewType,MatrixViewType,1>;
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
890  using FunctorType = FunctorRealSpaceTools::F_det<DeterminantArrayViewType,MatrixViewType,2>;
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
946  RealSpaceTools<DeviceType>::
947  add( Kokkos::DynRankView<sumArrayValueType,sumArrayProperties...> sumArray,
948  const Kokkos::DynRankView<inArray1ValueType,inArray1Properties...> inArray1,
949  const Kokkos::DynRankView<inArray2ValueType,inArray2Properties...> inArray2 ) {
950 
951 #ifdef HAVE_INTREPID2_DEBUG
952  {
953  INTREPID2_TEST_FOR_EXCEPTION( inArray1.rank() != inArray2.rank() ||
954  inArray1.rank() != sumArray.rank(), std::invalid_argument,
955  ">>> ERROR (RealSpaceTools::add): Array arguments must have identical ranks!");
956  for (size_type i=0;i<inArray1.rank();++i) {
957  INTREPID2_TEST_FOR_EXCEPTION( inArray1.extent(i) != inArray2.extent(i) ||
958  inArray1.extent(i) != sumArray.extent(i), std::invalid_argument,
959  ">>> ERROR (RealSpaceTools::add): Dimensions of array arguments do not agree!");
960  }
961  }
962 #endif
963  using MemSpaceType = typename DeviceType::memory_space;
964  constexpr bool are_accessible =
965  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
966  typename decltype(sumArray)::memory_space>::accessible &&
967  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
968  typename decltype(inArray1)::memory_space>::accessible &&
969  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
970  typename decltype(inArray2)::memory_space>::accessible;
971  static_assert(are_accessible, "RealSpaceTools<DeviceType>::add(..): input/output views' memory spaces are not compatible with DeviceType");
972 
973  using FunctorType = FunctorRealSpaceTools::F_add<decltype(sumArray),decltype(inArray1),decltype(inArray2)>;
974 
975  const auto loopSize = sumArray.extent(0);
976  using ExecSpaceType = typename DeviceType::execution_space;
977  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
978  Kokkos::parallel_for( policy, FunctorType(sumArray, inArray1, inArray2) );
979  }
980 
981  // ------------------------------------------------------------------------------------
982 
983  template<typename DeviceType>
984  template<typename inoutSumArrayValueType, class ...inoutSumArrayProperties,
985  typename inArrayValueType, class ...inArrayProperties>
986  void
987  RealSpaceTools<DeviceType>::
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
1045  RealSpaceTools<DeviceType>::
1046  subtract( Kokkos::DynRankView<diffArrayValueType,diffArrayProperties...> diffArray,
1047  const Kokkos::DynRankView<inArray1ValueType, inArray1Properties...> inArray1,
1048  const Kokkos::DynRankView<inArray2ValueType, inArray2Properties...> inArray2 ) {
1049 
1050 #ifdef HAVE_INTREPID2_DEBUG
1051  {
1052  INTREPID2_TEST_FOR_EXCEPTION( inArray1.rank() != inArray2.rank() ||
1053  inArray1.rank() != diffArray.rank(), std::invalid_argument,
1054  ">>> ERROR (RealSpaceTools::subtract): Array arguments must have identical ranks!");
1055  for (size_type i=0;i<inArray1.rank();++i) {
1056  INTREPID2_TEST_FOR_EXCEPTION( inArray1.extent(i) != inArray2.extent(i) ||
1057  inArray1.extent(i) != diffArray.extent(i), std::invalid_argument,
1058  ">>> ERROR (RealSpaceTools::subtract): Dimensions of array arguments do not agree!");
1059  }
1060  }
1061 #endif
1062  using MemSpaceType = typename DeviceType::memory_space;
1063  constexpr bool are_accessible =
1064  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1065  typename decltype(diffArray)::memory_space>::accessible &&
1066  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1067  typename decltype(inArray1)::memory_space>::accessible &&
1068  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1069  typename decltype(inArray2)::memory_space>::accessible;
1070  static_assert(are_accessible, "RealSpaceTools<DeviceType>::subtract(..): input/output views' memory spaces are not compatible with DeviceType");
1071 
1072  using FunctorType = FunctorRealSpaceTools::F_subtract<decltype(diffArray),decltype(inArray1),decltype(inArray2)>;
1073 
1074  const size_type loopSize = diffArray.extent(0);
1075  using ExecSpaceType = typename DeviceType::execution_space;
1076  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
1077  Kokkos::parallel_for( policy, FunctorType(diffArray, inArray1, inArray2) );
1078  }
1079 
1080  // ------------------------------------------------------------------------------------
1081 
1082  template<typename DeviceType>
1083  template<typename inoutDiffArrayValueType, class ...inoutDiffArrayProperties,
1084  typename inArrayValueType, class ...inArrayProperties>
1085  void
1086  RealSpaceTools<DeviceType>::
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
1144  RealSpaceTools<DeviceType>::
1145  scale( Kokkos::DynRankView<scaledArrayValueType,scaledArrayProperties...> scaledArray,
1146  const Kokkos::DynRankView<inArrayValueType, inArrayProperties...> inArray,
1147  const ValueType alpha ) {
1148 
1149 #ifdef HAVE_INTREPID2_DEBUG
1150  {
1151  INTREPID2_TEST_FOR_EXCEPTION( inArray.rank() > 5, std::invalid_argument,
1152  ">>> ERROR (RealSpaceTools::scale): Input array container has rank larger than 5.");
1153  INTREPID2_TEST_FOR_EXCEPTION( inArray.rank() != scaledArray.rank(), std::invalid_argument,
1154  ">>> ERROR (RealSpaceTools::scale): Array arguments must have identical ranks!");
1155  for (size_type i=0;i<inArray.rank();++i) {
1156  INTREPID2_TEST_FOR_EXCEPTION( inArray.extent(i) != scaledArray.extent(i), std::invalid_argument,
1157  ">>> ERROR (RealSpaceTools::scale): Dimensions of array arguments do not agree!");
1158  }
1159  }
1160 #endif
1161 
1162  using MemSpaceType = typename DeviceType::memory_space;
1163  constexpr bool are_accessible =
1164  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1165  typename decltype(scaledArray)::memory_space>::accessible &&
1166  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1167  typename decltype(inArray)::memory_space>::accessible;
1168  static_assert(are_accessible, "RealSpaceTools<DeviceType>::scale(..): input/output views' memory spaces are not compatible with DeviceType");
1169 
1170  using FunctorType = FunctorRealSpaceTools::F_scale<ValueType,decltype(scaledArray),decltype(inArray)>;
1171 
1172  const auto loopSize = scaledArray.extent(0);
1173  using ExecSpaceType = typename DeviceType::execution_space;
1174  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
1175  Kokkos::parallel_for( policy, FunctorType(scaledArray, inArray, alpha) );
1176  }
1177 
1178  // ------------------------------------------------------------------------------------
1179 
1180  template<typename DeviceType>
1181  template<typename ValueType,
1182  typename inoutScaledArrayValueType, class ...inoutScaledArrayProperties>
1183  void
1184  RealSpaceTools<DeviceType>::
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
1236  RealSpaceTools<DeviceType>::
1237  dot( Kokkos::DynRankView<dotArrayValueType,dotArrayProperties...> dotArray,
1238  const Kokkos::DynRankView<inVec1ValueType, inVec1Properties...> inVecs1,
1239  const Kokkos::DynRankView<inVec2ValueType, inVec2Properties...> inVecs2 ) {
1240 
1241 #ifdef HAVE_INTREPID2_DEBUG
1242  {
1243  INTREPID2_TEST_FOR_EXCEPTION( inVecs1.rank() != (dotArray.rank()+1), std::invalid_argument,
1244  ">>> ERROR (RealSpaceTools::dot): Ranks of norm and vector array arguments are incompatible!");
1245  INTREPID2_TEST_FOR_EXCEPTION( inVecs1.rank() != inVecs2.rank(), std::invalid_argument,
1246  ">>> ERROR (RealSpaceTools::dot): Ranks of input vector arguments must be identical!");
1247  INTREPID2_TEST_FOR_EXCEPTION( inVecs1.rank() < 2 || inVecs1.rank() > 3, std::invalid_argument,
1248  ">>> ERROR (RealSpaceTools::dot): Rank of input vector arguments must be 2 or 3!");
1249  for (size_type i=0;i<inVecs1.rank();++i) {
1250  INTREPID2_TEST_FOR_EXCEPTION( inVecs1.extent(i) != inVecs2.extent(i), std::invalid_argument,
1251  ">>> ERROR (RealSpaceTools::dot): Dimensions of input vector arguments do not agree!");
1252  }
1253  for (size_type i=0;i<dotArray.rank();++i) {
1254  INTREPID2_TEST_FOR_EXCEPTION( inVecs1.extent(i) != dotArray.extent(i), std::invalid_argument,
1255  ">>> ERROR (RealSpaceTools::dot): Dimensions of dot-product and vector arrays do not agree!");
1256  }
1257  }
1258 #endif
1259  using MemSpaceType = typename DeviceType::memory_space;
1260  constexpr bool are_accessible =
1261  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1262  typename decltype(dotArray)::memory_space>::accessible &&
1263  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1264  typename decltype(inVecs1)::memory_space>::accessible &&
1265  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1266  typename decltype(inVecs2)::memory_space>::accessible;
1267  static_assert(are_accessible, "RealSpaceTools<DeviceType>::dot(..): input/output views' memory spaces are not compatible with DeviceType");
1268 
1269  using FunctorType = FunctorRealSpaceTools::F_dot<decltype(dotArray),decltype(inVecs1),decltype(inVecs2)>;
1270 
1271  const auto loopSize = dotArray.extent(0)*dotArray.extent(1);
1272  using ExecSpaceType = typename DeviceType::execution_space;
1273  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
1274  Kokkos::parallel_for( policy, FunctorType(dotArray, inVecs1, inVecs2) );
1275  }
1276 
1277  // ------------------------------------------------------------------------------------
1278 
1279  namespace FunctorRealSpaceTools {
1283  template<typename matVecViewType,
1284  typename inMatViewType,
1285  typename inVecViewType>
1286  struct F_matvec {
1287  matVecViewType _matVecs;
1288  const inMatViewType _inMats;
1289  const inVecViewType _inVecs;
1290 
1291  KOKKOS_INLINE_FUNCTION
1292  F_matvec( matVecViewType matVecs_,
1293  inMatViewType inMats_,
1294  inVecViewType inVecs_ )
1295  : _matVecs(matVecs_), _inMats(inMats_), _inVecs(inVecs_) {}
1296 
1297  KOKKOS_INLINE_FUNCTION
1298  void operator()(const ordinal_type iter) const {
1299  const ordinal_type rm = _inMats.rank(), rv = _inVecs.rank(), rr = _matVecs.rank();
1300  ordinal_type _i = iter, _j = 0;
1301 
1302  if ( rr > 2 )
1303  unrollIndex( _i, _j,
1304  _matVecs.extent(0),
1305  _matVecs.extent(1),
1306  iter );
1307 
1308  auto mat = ( rm == 2 ? Kokkos::subview(_inMats, Kokkos::ALL(), Kokkos::ALL()) :
1309  rm == 3 ? Kokkos::subview(_inMats, _i, Kokkos::ALL(), Kokkos::ALL()) :
1310  Kokkos::subview(_inMats, _i, _j, Kokkos::ALL(), Kokkos::ALL()) );
1311 
1312  auto vec = ( rv == 1 ? Kokkos::subview(_inVecs, Kokkos::ALL()) :
1313  rv == 2 ? Kokkos::subview(_inVecs, _i, Kokkos::ALL()) :
1314  Kokkos::subview(_inVecs, _i, _j, Kokkos::ALL()) );
1315 
1316  auto result = ( rr == 1 ? Kokkos::subview(_matVecs, Kokkos::ALL()) :
1317  rr == 2 ? Kokkos::subview(_matVecs, _i, Kokkos::ALL()) :
1318  Kokkos::subview(_matVecs, _i, _j, Kokkos::ALL()) );
1319 
1320  const ordinal_type iend = result.extent(0);
1321  const ordinal_type jend = vec.extent(0);
1322 
1323  for (ordinal_type i=0;i<iend;++i) {
1324  result(i) = 0;
1325  for (ordinal_type j=0;j<jend;++j)
1326  result(i) += mat(i, j)*vec(j);
1327  }
1328  }
1329  };
1330 
1331 
1335  template<typename outMatViewType,
1336  typename inMatViewType>
1337  struct F_AtA {
1338  outMatViewType _outMats;
1339  const inMatViewType _inMats;
1340 
1341  KOKKOS_INLINE_FUNCTION
1342  F_AtA( outMatViewType outMats_,
1343  inMatViewType inMats_)
1344  : _outMats(outMats_), _inMats(inMats_) {}
1345 
1346  KOKKOS_INLINE_FUNCTION
1347  void operator()(const ordinal_type iter) const {
1348  const ordinal_type rIM = _inMats.rank(), rOM = _outMats.rank();
1349  ordinal_type _i = iter, _j = 0;
1350 
1351  if ( rOM > 3 )
1352  unrollIndex( _i, _j,
1353  _outMats.extent(0),
1354  _outMats.extent(1),
1355  iter );
1356 
1357  auto inMat = ( rIM == 2 ? Kokkos::subview(_inMats, Kokkos::ALL(), Kokkos::ALL()) :
1358  rIM == 3 ? Kokkos::subview(_inMats, _i, Kokkos::ALL(), Kokkos::ALL()) :
1359  Kokkos::subview(_inMats, _i, _j, Kokkos::ALL(), Kokkos::ALL()) );
1360 
1361  auto outMat = ( rOM == 2 ? Kokkos::subview(_outMats, Kokkos::ALL(), Kokkos::ALL()) :
1362  rOM == 3 ? Kokkos::subview(_outMats, _i, Kokkos::ALL(), Kokkos::ALL()) :
1363  Kokkos::subview(_outMats, _i, _j, Kokkos::ALL(), Kokkos::ALL()) );
1364 
1365  const ordinal_type iend = outMat.extent(0);
1366  const ordinal_type jend = outMat.extent(1);
1367  const ordinal_type kend = inMat.extent(0);
1368 
1369  for (ordinal_type i=0;i<iend;++i) {
1370  for (ordinal_type j=i;j<jend;++j) {
1371  outMat(i,j) = 0;
1372  for (ordinal_type k=0;k<kend;++k)
1373  outMat(i,j) += inMat(k, i)*inMat(k,j);
1374  }
1375  for (ordinal_type j=0;j<i;++j)
1376  outMat(i,j) = outMat(j,i);
1377  }
1378  }
1379  };
1380 
1381  }
1382 
1383  template<typename DeviceType>
1384  template<typename matVecValueType, class ...matVecProperties,
1385  typename inMatValueType, class ...inMatProperties,
1386  typename inVecValueType, class ...inVecProperties>
1387  void
1388  RealSpaceTools<DeviceType>::
1389  matvec( Kokkos::DynRankView<matVecValueType,matVecProperties...> matVecs,
1390  const Kokkos::DynRankView<inMatValueType, inMatProperties...> inMats,
1391  const Kokkos::DynRankView<inVecValueType, inVecProperties...> inVecs ) {
1392 
1393 #ifdef HAVE_INTREPID2_DEBUG
1394  INTREPID2_TEST_FOR_EXCEPTION( inMats.rank() < 2 || inMats.rank() > 4, std::invalid_argument,
1395  ">>> ERROR (RealSpaceTools::matvec): Rank of matrix array must be 2, 3 or 4!");
1396  INTREPID2_TEST_FOR_EXCEPTION( matVecs.rank() < 1 || matVecs.rank() > 3, std::invalid_argument,
1397  ">>> ERROR (RealSpaceTools::matvec): Rank of matVecs array must be 1, 2 or 3!");
1398  INTREPID2_TEST_FOR_EXCEPTION( inVecs.rank() < 1 || inVecs.rank() > 3, std::invalid_argument,
1399  ">>> ERROR (RealSpaceTools::matvec): Rank of inVecs array must be 1, 2 or 3!");
1400  if (inMats.rank() == 2) {
1401  // a single matrix, multiple input and output
1402  INTREPID2_TEST_FOR_EXCEPTION( matVecs.rank() != inVecs.rank(), std::invalid_argument,
1403  ">>> ERROR (RealSpaceTools::matvec): Vector arrays must be have the same rank!");
1404  // output must match
1405  for (ordinal_type i=0;i< (static_cast<ordinal_type>(inVecs.rank())-1);++i) {
1406  INTREPID2_TEST_FOR_EXCEPTION( matVecs.extent(i) != inVecs.extent(i), std::invalid_argument,
1407  ">>> ERROR (RealSpaceTools::matvec): Dimensions of vector array arguments do not agree!");
1408  }
1409  // matvec compatibility
1410  INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(0) != matVecs.extent(matVecs.rank()-1) ||
1411  inMats.extent(1) != inVecs.extent(inVecs.rank()-1), std::invalid_argument,
1412  ">>> ERROR (RealSpaceTools::matvec): Matvec dimensions are not compatible each other.");
1413  } else if (inVecs.rank() < matVecs.rank()) {
1414  // multiple matrix, single input and multiple output
1415  INTREPID2_TEST_FOR_EXCEPTION( inMats.rank() != (matVecs.rank()+1), std::invalid_argument,
1416  ">>> ERROR (RealSpaceTools::matvec): The result vector and matrix array arguments do not have compatible ranks!");
1417  for (ordinal_type i=0;i<(static_cast<ordinal_type>(inMats.rank())-2);++i) {
1418  INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(i) != matVecs.extent(i), std::invalid_argument,
1419  ">>> ERROR (RealSpaceTools::matvec): Dimensions of vector and matrix array arguments do not agree!");
1420  }
1421  // matvec compatibility
1422  INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(inMats.rank()-2) != matVecs.extent(matVecs.rank()-1) ||
1423  inMats.extent(inMats.rank()-1) != inVecs.extent(inVecs.rank()-1), std::invalid_argument,
1424  ">>> ERROR (RealSpaceTools::matvec): Matvec dimensions are not compatible each other.");
1425  } else {
1426  // multiple matrix, multiple input and multiple output
1427  INTREPID2_TEST_FOR_EXCEPTION( inMats.rank() != (matVecs.rank()+1), std::invalid_argument,
1428  ">>> ERROR (RealSpaceTools::matvec): The result vector and matrix array arguments do not have compatible ranks!");
1429  INTREPID2_TEST_FOR_EXCEPTION( matVecs.rank() != inVecs.rank(), std::invalid_argument,
1430  ">>> ERROR (RealSpaceTools::matvec): Vector arrays must be have the same rank!");
1431  for (ordinal_type i=0;i<(static_cast<ordinal_type>(inMats.rank())-2);++i) {
1432  INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(i) != matVecs.extent(i), std::invalid_argument,
1433  ">>> ERROR (RealSpaceTools::matvec): Dimensions of vector and matrix array arguments do not agree!");
1434  }
1435  for (size_type i=0;i<inVecs.rank();++i) {
1436  INTREPID2_TEST_FOR_EXCEPTION( matVecs.extent(i) != inVecs.extent(i), std::invalid_argument,
1437  ">>> ERROR (RealSpaceTools::matvec): Dimensions of vector array arguments do not agree!");
1438  }
1439  INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(inMats.rank()-1) != inVecs.extent(inVecs.rank()-1), std::invalid_argument,
1440  ">>> ERROR (RealSpaceTools::matvec): Matrix column dimension does not match to the length of a vector!");
1441  }
1442 #endif
1443  using MemSpaceType = typename DeviceType::memory_space;
1444  constexpr bool are_accessible =
1445  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1446  typename decltype(matVecs)::memory_space>::accessible &&
1447  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1448  typename decltype(inMats)::memory_space>::accessible &&
1449  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1450  typename decltype(inVecs)::memory_space>::accessible;
1451  static_assert(are_accessible, "RealSpaceTools<DeviceType>::matvec(..): input/output views' memory spaces are not compatible with DeviceType");
1452 
1453  using FunctorType = FunctorRealSpaceTools::F_matvec<decltype(matVecs),decltype(inMats),decltype(inVecs)>;
1454 
1455  size_type loopSize = 1;
1456  const ordinal_type r = matVecs.rank() - 1;
1457  for (ordinal_type i=0;i<r;++i)
1458  loopSize *= matVecs.extent(i);
1459 
1460  using ExecSpaceType = typename DeviceType::execution_space;
1461  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
1462  Kokkos::parallel_for( policy, FunctorType(matVecs, inMats, inVecs) );
1463  }
1464 
1465 
1466  template<typename DeviceType>
1467  template<typename outMatValueType, class ...outMatProperties,
1468  typename inMatValueType, class ...inMatProperties>
1469  void
1470  RealSpaceTools<DeviceType>::
1471  AtA( Kokkos::DynRankView<outMatValueType,outMatProperties...> outMats,
1472  const Kokkos::DynRankView<inMatValueType, inMatProperties...> inMats) {
1473 
1474 #ifdef HAVE_INTREPID2_DEBUG
1475  INTREPID2_TEST_FOR_EXCEPTION( inMats.rank() < 2 || inMats.rank() > 4, std::invalid_argument,
1476  ">>> ERROR (RealSpaceTools::AtA): Rank of input matrix array must be 2, 3 or 4!");
1477 
1478  INTREPID2_TEST_FOR_EXCEPTION( inMats.rank() != outMats.rank(), std::invalid_argument,
1479  ">>> ERROR (RealSpaceTools::AtA): The matrices do not have the same ranks!");
1480  for (ordinal_type i=0;i<(static_cast<ordinal_type>(inMats.rank())-2);++i) {
1481  INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(i) != outMats.extent(i), std::invalid_argument,
1482  ">>> ERROR (RealSpaceTools::AtA): Dimensions of matrices arrays do not agree!");
1483  }
1484  // matmat compatibility
1485  INTREPID2_TEST_FOR_EXCEPTION( inMats.extent(inMats.rank()-1) != outMats.extent(outMats.rank()-1) ||
1486  outMats.extent(outMats.rank()-1) != outMats.extent(outMats.rank()-2), std::invalid_argument,
1487  ">>> ERROR (RealSpaceTools::AtA): Matrices dimensions are not compatible.");
1488 
1489 #endif
1490  using MemSpaceType = typename DeviceType::memory_space;
1491  constexpr bool are_accessible =
1492  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1493  typename decltype(outMats)::memory_space>::accessible &&
1494  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1495  typename decltype(inMats)::memory_space>::accessible;
1496  static_assert(are_accessible, "RealSpaceTools<DeviceType>::AtA(..): input/output views' memory spaces are not compatible with DeviceType");
1497 
1498  using FunctorType = FunctorRealSpaceTools::F_AtA<decltype(outMats),decltype(inMats)>;
1499 
1500  size_type loopSize = 1;
1501  const ordinal_type r = outMats.rank() - 2;
1502  for (ordinal_type i=0;i<r;++i)
1503  loopSize *= outMats.extent(i);
1504 
1505  using ExecSpaceType = typename DeviceType::execution_space;
1506  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
1507  Kokkos::parallel_for( policy, FunctorType(outMats, inMats) );
1508  }
1509 
1510 
1511  // ------------------------------------------------------------------------------------
1512 
1513  namespace FunctorRealSpaceTools {
1517  template<typename vecProdViewType,
1518  typename inLeftViewType,
1519  typename inRightViewType>
1520  struct F_vecprod {
1521  vecProdViewType _vecProd;
1522  const inLeftViewType _inLeft;
1523  const inRightViewType _inRight;
1524  const bool _is_vecprod_3d;
1525 
1526  KOKKOS_INLINE_FUNCTION
1527  F_vecprod( vecProdViewType vecProd_,
1528  inLeftViewType inLeft_,
1529  inRightViewType inRight_,
1530  const bool is_vecprod_3d_ )
1531  : _vecProd(vecProd_), _inLeft(inLeft_), _inRight(inRight_), _is_vecprod_3d(is_vecprod_3d_) {}
1532 
1533  KOKKOS_INLINE_FUNCTION
1534  void operator()(const size_type iter) const {
1535  ordinal_type i, j;
1536  unrollIndex( i, j,
1537  _inLeft.extent(0),
1538  _inLeft.extent(1),
1539  iter );
1540 
1541  const auto r = _inLeft.rank();
1542 
1543  auto left = ( r == 1 ? Kokkos::subview(_inLeft, Kokkos::ALL()) :
1544  r == 2 ? Kokkos::subview(_inLeft, i, Kokkos::ALL()) :
1545  Kokkos::subview(_inLeft, i, j, Kokkos::ALL()) );
1546 
1547  auto right = ( r == 1 ? Kokkos::subview(_inRight, Kokkos::ALL()) :
1548  r == 2 ? Kokkos::subview(_inRight, i, Kokkos::ALL()) :
1549  Kokkos::subview(_inRight, i, j, Kokkos::ALL()) );
1550 
1551  auto result = ( r == 1 ? Kokkos::subview(_vecProd, Kokkos::ALL()) :
1552  r == 2 ? Kokkos::subview(_vecProd, i, Kokkos::ALL()) :
1553  Kokkos::subview(_vecProd, i, j, Kokkos::ALL()) );
1554 
1555  // branch prediction is possible
1556  if (_is_vecprod_3d) {
1557  result(0) = left(1)*right(2) - left(2)*right(1);
1558  result(1) = left(2)*right(0) - left(0)*right(2);
1559  result(2) = left(0)*right(1) - left(1)*right(0);
1560  } else {
1561  result(0) = left(0)*right(1) - left(1)*right(0);
1562  }
1563  }
1564  };
1565  }
1566 
1567  template<typename DeviceType>
1568  template<typename vecProdValueType, class ...vecProdProperties,
1569  typename inLeftValueType, class ...inLeftProperties,
1570  typename inRightValueType, class ...inRightProperties>
1571  void
1572  RealSpaceTools<DeviceType>::
1573  vecprod( Kokkos::DynRankView<vecProdValueType,vecProdProperties...> vecProd,
1574  const Kokkos::DynRankView<inLeftValueType, inLeftProperties...> inLeft,
1575  const Kokkos::DynRankView<inRightValueType,inRightProperties...> inRight ) {
1576 
1577 #ifdef HAVE_INTREPID2_DEBUG
1578  {
1579  /*
1580  * Check array rank and spatial dimension range (if applicable)
1581  * (1) all array arguments are required to have matching dimensions and rank: (D), (I0,D) or (I0,I1,D)
1582  * (2) spatial dimension should be 2 or 3
1583  */
1584  // (1) check rank range on inLeft and then compare the other arrays with inLeft
1585  INTREPID2_TEST_FOR_EXCEPTION( inLeft.rank() < 1 || inLeft.rank() > 3, std::invalid_argument,
1586  ">>> ERROR (RealSpaceTools::vecprod): Rank of matrix array must be 1, 2, or 3!");
1587  INTREPID2_TEST_FOR_EXCEPTION( inLeft.rank() != inRight.rank(), std::invalid_argument,
1588  ">>> ERROR (RealSpaceTools::vecprod): Right and left arrays must be have the same rank!");
1589  INTREPID2_TEST_FOR_EXCEPTION( inLeft.rank() != vecProd.rank(), std::invalid_argument,
1590  ">>> ERROR (RealSpaceTools::vecprod): Left and vecProd arrays must be have the same rank!");
1591  for (size_type i=0;i<inLeft.rank();++i) {
1592  INTREPID2_TEST_FOR_EXCEPTION( inLeft.extent(i) != inRight.extent(i), std::invalid_argument,
1593  ">>> ERROR (RealSpaceTools::vecprod): Dimensions of matrix arguments do not agree!");
1594  INTREPID2_TEST_FOR_EXCEPTION( inLeft.extent(i) != vecProd.extent(i), std::invalid_argument,
1595  ">>> ERROR (RealSpaceTools::vecprod): Dimensions of matrix arguments do not agree!");
1596  }
1597 
1598  // (2) spatial dimension ordinal = array rank - 1. Suffices to check only one array because we just
1599  // checked whether or not the arrays have matching dimensions.
1600  INTREPID2_TEST_FOR_EXCEPTION( inLeft.extent(inLeft.rank()-1) < 2 ||
1601  inLeft.extent(inLeft.rank()-1) > 3, std::invalid_argument,
1602  ">>> ERROR (RealSpaceTools::vecprod): Dimensions of arrays (rank-1) must be 2 or 3!");
1603  }
1604 #endif
1605  using MemSpaceType = typename DeviceType::memory_space;
1606  constexpr bool are_accessible =
1607  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1608  typename decltype(vecProd)::memory_space>::accessible &&
1609  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1610  typename decltype(inLeft)::memory_space>::accessible &&
1611  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
1612  typename decltype(inRight)::memory_space>::accessible;
1613  static_assert(are_accessible, "RealSpaceTools<DeviceType>::vecprod(..): input/output views' memory spaces are not compatible with DeviceType");
1614 
1615  using FunctorType = FunctorRealSpaceTools::F_vecprod<decltype(vecProd),decltype(inLeft),decltype(inRight)>;
1616 
1617  const auto r = inLeft.rank();
1618  const auto loopSize = ( r == 1 ? 1 :
1619  r == 2 ? inLeft.extent(0) :
1620  inLeft.extent(0)*inLeft.extent(1) );
1621  const bool is_vecprod_3d = (inLeft.extent(inLeft.rank() - 1) == 3);
1622  using ExecSpaceType = typename DeviceType::execution_space;
1623  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
1624  Kokkos::parallel_for( policy, FunctorType(vecProd, inLeft, inRight, is_vecprod_3d) );
1625  }
1626 
1627  // ------------------------------------------------------------------------------------
1628 
1629 } // namespace Intrepid2
1630 
1631 #endif
1632 
1633 
1634 
1635 
1636 
1637 // =====================================
1638 // Too much things...
1639 //
1640 
1641 // template<class ...inVec1Properties,
1642 // class ...inVec2Properties>
1643 // KOKKOS_INLINE_FUNCTION
1644 // typename Kokkos::DynRankView<inVec1Properties...>::value_type
1645 // RealSpaceTools<DeviceType>::
1646 // dot( const Kokkos::DynRankView<inVec1Properties...> inVec1,
1647 // const Kokkos::DynRankView<inVec2Properties...> inVec2 ) {
1648 
1649 // #ifdef HAVE_INTREPID2_DEBUG
1650 // INTREPID2_TEST_FOR_ABORT( inVec1.rank != 1 || inVec2.rank() != 1,
1651 // ">>> ERROR (RealSpaceTools::dot): Vector arguments must have rank 1!");
1652 // INTREPID2_TEST_FOR_ABORT( inVec1.extent(0) != inVec2.extent(0),
1653 // ">>> ERROR (RealSpaceTools::dot): Dimensions of vector arguments must agree!");
1654 // #endif
1655 // typedef typename Kokkos::DynRankView<inVec1Properties...>::value_type value_type;
1656 
1657 // // designed for small problems
1658 // value_type r_val(0);
1659 
1660 // // * This is probably not necessary
1661 // if (std::is_same<ExecSpace,Kokkos::Serial>::value) {
1662 // const ordinal_type iend = iVec1.extent(0);
1663 // for (ordinal_type i=0;i<iend;++i)
1664 // r_val += inVec1(i)*inVec2(i);
1665 // } else {
1666 // struct Functor {
1667 // Kokkos::DynRankView<inVec1Properties...> _inVec1;
1668 // Kokkos::DynRankView<inVec2Properties...> _inVec2;
1669 
1670 // KOKKOS_INLINE_FUNCTION
1671 // Functor(const Kokkos::DynRankView<inVec1Properties...> inVec1_,
1672 // const Kokkos::DynRankView<inVec2Properties...> inVec2_);
1673 
1674 // KOKKOS_INLINE_FUNCTION
1675 // void operator()(ordinal_type i, value_type &dst ) const {
1676 // dst += _inVec1(i)*_inVec2(i);
1677 // }
1678 
1679 // KOKKOS_INLINE_FUNCTION
1680 // void join(volatile value_type &dst ,
1681 // const volatile value_type &src) const {
1682 // dst += src;
1683 // }
1684 // };
1685 // const ordinal_type iend = inVec1.extent(0);
1686 // Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, iend);
1687 // Kokkos::parallel_for( policy, Functor(inVec1, inVec2), r_val );
1688 // }
1689 // return r_val;
1690 // }
1691 
1692 
Intrepid2::RealSpaceTools::transpose
static void transpose(Kokkos::DynRankView< transposeMatValueType, transposeMatProperties...> transposeMats, const Kokkos::DynRankView< inMatValueType, inMatProperties...> inMats)
Computes transposes of square matrices stored in an array of total rank 2 (single matrix)...
Definition: Intrepid2_RealSpaceToolsDef.hpp:585
Intrepid2::RealSpaceTools::vectorNorm
static void vectorNorm(Kokkos::DynRankView< normArrayValueType, normArrayProperties...> normArray, const Kokkos::DynRankView< inVecValueType, inVecProperties...> inVecs, const ENorm normType)
Computes norms (1, 2, infinity) of vectors stored in a array of total rank 2 (array of vectors)...
Definition: Intrepid2_RealSpaceToolsDef.hpp:500
Intrepid2::RealSpaceTools::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:856
Intrepid2::Util
small utility functions
Definition: Intrepid2_Utils.hpp:270
Intrepid2::FunctorRealSpaceTools::F_dot
Functor to compute dot product see Intrepid2::RealSpaceTools for more.
Definition: Intrepid2_RealSpaceToolsDef.hpp:1200
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:1471
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:67
Intrepid2::RealSpaceTools::extractScalarValues
static void extractScalarValues(Kokkos::DynRankView< outputValueType, outputProperties...> output, const Kokkos::DynRankView< inputValueType, inputProperties...> input)
Extract scalar type values from Sacado-based array.
Definition: Intrepid2_RealSpaceToolsDef.hpp:257
Intrepid2::RealSpaceTools::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:1573
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:1046
Intrepid2::FunctorRealSpaceTools::F_inverse
Functor to compute inverse see Intrepid2::RealSpaceTools for more.
Definition: Intrepid2_RealSpaceToolsDef.hpp:631
Intrepid2::FunctorRealSpaceTools::F_det
Functor to compute determinant see Intrepid2::RealSpaceTools for more.
Definition: Intrepid2_RealSpaceToolsDef.hpp:802
Intrepid2::RealSpaceTools::clone
static void clone(Kokkos::DynRankView< outputValueType, outputProperties...> output, const Kokkos::DynRankView< inputValueType, inputProperties...> input)
Clone input array.
Definition: Intrepid2_RealSpaceToolsDef.hpp:334
Intrepid2::RealSpaceTools::Serial::dot
static KOKKOS_INLINE_FUNCTION inVec1ValueType dot(const Kokkos::DynRankView< inVec1ValueType, inVec1Properties...> inVec1, const Kokkos::DynRankView< inVec2ValueType, inVec2Properties...> inVec2)
Computes dot product of two vectors stored in arrays of rank 1.
Definition: Intrepid2_RealSpaceToolsDef.hpp:185
Intrepid2::RealSpaceTools::inverse
static void inverse(InverseMatrixViewType inverseMats, MatrixViewType inMats)
Computes inverses of nonsingular matrices stored in an array of total rank 2 (single matrix)...
Definition: Intrepid2_RealSpaceToolsDef.hpp:750
Intrepid2::FunctorRealSpaceTools::F_vecprod
Functor to compute vecprod see Intrepid2::RealSpaceTools for more.
Definition: Intrepid2_RealSpaceToolsDef.hpp:1520
Intrepid2::FunctorRealSpaceTools::F_absval
Functor to compute absolute value see Intrepid2::RealSpaceTools for more.
Definition: Intrepid2_RealSpaceToolsDef.hpp:387
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:1145
Intrepid2::FunctorRealSpaceTools::F_subtract
Functor to subtract md arrays see Intrepid2::RealSpaceTools for more.
Definition: Intrepid2_RealSpaceToolsDef.hpp:1013
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:137
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:1237
Intrepid2::FunctorRealSpaceTools::F_scale
Functor to scale md arrays see Intrepid2::RealSpaceTools for more.
Definition: Intrepid2_RealSpaceToolsDef.hpp:1111
Intrepid2::FunctorRealSpaceTools::F_add
Functor to add md arrays see Intrepid2::RealSpaceTools for more.
Definition: Intrepid2_RealSpaceToolsDef.hpp:914
Intrepid2::FunctorRealSpaceTools::F_vectorNorm
Functor to compute vector norm see Intrepid2::RealSpaceTools for more.
Definition: Intrepid2_RealSpaceToolsDef.hpp:468
Intrepid2::FunctorRealSpaceTools::F_extractScalarValues
Functor for extractScalarValues see Intrepid2::RealSpaceTools for more.
Definition: Intrepid2_RealSpaceToolsDef.hpp:227
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:1389
Intrepid2::FunctorRealSpaceTools::F_transpose
Functor to compute transpose see Intrepid2::RealSpaceTools for more.
Definition: Intrepid2_RealSpaceToolsDef.hpp:540
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:947
Intrepid2::FunctorRealSpaceTools::F_clone
Functor for clone see Intrepid2::RealSpaceTools for more.
Definition: Intrepid2_RealSpaceToolsDef.hpp:282
Intrepid2::RealSpaceTools::absval
static void absval(Kokkos::DynRankView< absArrayValueType, absArrayProperties...> absArray, const Kokkos::DynRankView< inArrayValueType, inArrayProperties...> inArray)
Computes absolute value of an array.
Definition: Intrepid2_RealSpaceToolsDef.hpp:419
Intrepid2::FunctorRealSpaceTools::F_matvec
Functor to compute matvec see Intrepid2::RealSpaceTools for more.
Definition: Intrepid2_RealSpaceToolsDef.hpp:1286
Intrepid2::FunctorRealSpaceTools::F_AtA
Functor to compute matvec see Intrepid2::RealSpaceTools for more.
Definition: Intrepid2_RealSpaceToolsDef.hpp:1337
Generated by   doxygen 1.8.5