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