Intrepid
test_01_kokkos.cpp
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 Pavel Bochev (pbboche@sandia.gov)
38 // Denis Ridzal (dridzal@sandia.gov), or
39 // Kara Peterson (kjpeter@sandia.gov)
40 //
41 // ************************************************************************
42 // @HEADER
43 
49 #include "Intrepid_ArrayTools.hpp"
50 #include "Intrepid_FieldContainer.hpp"
51 #include "Intrepid_RealSpaceTools.hpp"
52 #include "Teuchos_oblackholestream.hpp"
53 #include "Teuchos_RCP.hpp"
54 #include "Teuchos_ScalarTraits.hpp"
55 #include "Teuchos_GlobalMPISession.hpp"
56 #include <Kokkos_Core.hpp>
57 
58 using namespace std;
59 using namespace Intrepid;
60 
61 #define INTREPID_TEST_COMMAND( S ) \
62 { \
63  try { \
64  S ; \
65  } \
66  catch (std::logic_error err) { \
67  *outStream << "Expected Error ----------------------------------------------------------------\n"; \
68  *outStream << err.what() << '\n'; \
69  *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
70  }; \
71 }
72 
73 
74 int main(int argc, char *argv[]) {
75 
76  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
77  Kokkos::initialize();
78  // This little trick lets us print to std::cout only if
79  // a (dummy) command-line argument is provided.
80  int iprint = argc - 1;
81  Teuchos::RCP<std::ostream> outStream;
82  Teuchos::oblackholestream bhs; // outputs nothing
83  if (iprint > 0)
84  outStream = Teuchos::rcp(&std::cout, false);
85  else
86  outStream = Teuchos::rcp(&bhs, false);
87 
88  // Save the format state of the original std::cout.
89  Teuchos::oblackholestream oldFormatState;
90  oldFormatState.copyfmt(std::cout);
91 
92  *outStream \
93  << "===============================================================================\n" \
94  << "| |\n" \
95  << "| Unit Test (ArrayTools) |\n" \
96  << "| |\n" \
97  << "| 1) Array operations: contractions |\n" \
98  << "| |\n" \
99  << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \
100  << "| Denis Ridzal (dridzal@sandia.gov). |\n" \
101  << "| |\n" \
102  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
103  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
104  << "| |\n" \
105  << "===============================================================================\n";
106 
107 
108  int errorFlag = 0;
109 
110 #ifdef HAVE_INTREPID_DEBUG
111  int beginThrowNumber = Teuchos::TestForException_getThrowNumber();
112  int endThrowNumber = beginThrowNumber + 81;
113 #endif
114  typedef ArrayTools art;
115  typedef RealSpaceTools<double> rst;
116 #ifdef HAVE_INTREPID_DEBUG
117  ArrayTools atools;
118 #endif
119  *outStream \
120  << "\n"
121  << "===============================================================================\n"\
122  << "| TEST 1: exceptions |\n"\
123  << "===============================================================================\n";
124 
125  try{
126 
127 #ifdef HAVE_INTREPID_DEBUG
128  FieldContainer<double> a_2_2(2, 2);
129  FieldContainer<double> a_10_2(10, 2);
130  FieldContainer<double> a_10_3(10, 3);
131  FieldContainer<double> a_10_2_2(10, 2, 2);
132  FieldContainer<double> a_10_2_3(10, 2, 3);
133  FieldContainer<double> a_10_3_2(10, 3, 2);
134  FieldContainer<double> a_9_2_2(9, 2, 2);
135 
136  *outStream << "-> contractFieldFieldScalar:\n";
137  INTREPID_TEST_COMMAND( atools.contractFieldFieldScalar<double>(a_2_2, a_2_2, a_2_2, COMP_CPP) );
138  INTREPID_TEST_COMMAND( atools.contractFieldFieldScalar<double>(a_2_2, a_10_2_2, a_2_2, COMP_CPP) );
139  INTREPID_TEST_COMMAND( atools.contractFieldFieldScalar<double>(a_2_2, a_10_2_2, a_10_2_2, COMP_CPP) );
140  INTREPID_TEST_COMMAND( atools.contractFieldFieldScalar<double>(a_10_2_2, a_9_2_2, a_10_2_2, COMP_CPP) );
141  INTREPID_TEST_COMMAND( atools.contractFieldFieldScalar<double>(a_10_2_2, a_10_2_2, a_10_2_3, COMP_CPP) );
142  INTREPID_TEST_COMMAND( atools.contractFieldFieldScalar<double>(a_9_2_2, a_10_3_2, a_10_2_2, COMP_CPP) );
143  INTREPID_TEST_COMMAND( atools.contractFieldFieldScalar<double>(a_10_3_2, a_10_2_2, a_10_3_2, COMP_CPP) );
144  INTREPID_TEST_COMMAND( atools.contractFieldFieldScalar<double>(a_10_2_2, a_10_2_2, a_10_3_2, COMP_CPP) );
145  INTREPID_TEST_COMMAND( atools.contractFieldFieldScalar<double>(a_10_2_3, a_10_2_2, a_10_3_2, COMP_ENGINE_MAX) );
146  INTREPID_TEST_COMMAND( atools.contractFieldFieldScalar<double>(a_10_2_3, a_10_2_2, a_10_3_2, COMP_CPP) );
147 
148 
149  FieldContainer<double> a_10_2_2_2(10, 2, 2, 2);
150  FieldContainer<double> a_9_2_2_2(9, 2, 2, 2);
151  FieldContainer<double> a_10_3_2_2(10, 3, 2, 2);
152  FieldContainer<double> a_10_2_3_2(10, 2, 3, 2);
153  FieldContainer<double> a_10_2_2_3(10, 2, 2, 3);
154 
155 
156  *outStream << "-> contractFieldFieldVector:\n";
157  INTREPID_TEST_COMMAND( atools.contractFieldFieldVector<double>(a_2_2, a_2_2, a_2_2, COMP_CPP) );
158  INTREPID_TEST_COMMAND( atools.contractFieldFieldVector<double>(a_2_2, a_10_2_2_2, a_2_2, COMP_CPP) );
159  INTREPID_TEST_COMMAND( atools.contractFieldFieldVector<double>(a_2_2, a_10_2_2_2, a_10_2_2_2, COMP_CPP) );
160  INTREPID_TEST_COMMAND( atools.contractFieldFieldVector<double>(a_10_2_2, a_9_2_2_2, a_10_2_2_2, COMP_CPP) );
161  INTREPID_TEST_COMMAND( atools.contractFieldFieldVector<double>(a_10_2_2, a_10_2_2_2, a_10_2_3_2, COMP_CPP) );
162  INTREPID_TEST_COMMAND( atools.contractFieldFieldVector<double>(a_10_2_2, a_10_2_2_2, a_10_2_2_3, COMP_CPP) );
163  INTREPID_TEST_COMMAND( atools.contractFieldFieldVector<double>(a_9_2_2, a_10_2_2_2, a_10_2_2_2, COMP_CPP) );
164  INTREPID_TEST_COMMAND( atools.contractFieldFieldVector<double>(a_10_3_2, a_10_2_2_2, a_10_3_2_2, COMP_CPP) );
165  INTREPID_TEST_COMMAND( atools.contractFieldFieldVector<double>(a_10_2_2, a_10_2_2_2, a_10_3_2_2, COMP_CPP) );
166  INTREPID_TEST_COMMAND( atools.contractFieldFieldVector<double>(a_10_2_3, a_10_2_2_2, a_10_3_2_2, COMP_ENGINE_MAX) );
167  INTREPID_TEST_COMMAND( atools.contractFieldFieldVector<double>(a_10_2_3, a_10_2_2_2, a_10_3_2_2, COMP_CPP) );
168 
169 
170  FieldContainer<double> a_10_2_2_2_2(10, 2, 2, 2, 2);
171  FieldContainer<double> a_9_2_2_2_2(9, 2, 2, 2, 2);
172  FieldContainer<double> a_10_3_2_2_2(10, 3, 2, 2, 2);
173  FieldContainer<double> a_10_2_3_2_2(10, 2, 3, 2, 2);
174  FieldContainer<double> a_10_2_2_3_2(10, 2, 2, 3, 2);
175  FieldContainer<double> a_10_2_2_2_3(10, 2, 2, 2, 3);
176 
177  *outStream << "-> contractFieldFieldTensor:\n";
178  INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<double>(a_2_2, a_2_2, a_2_2, COMP_CPP) );
179  INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<double>(a_2_2, a_10_2_2_2_2, a_2_2, COMP_CPP) );
180  INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<double>(a_2_2, a_10_2_2_2_2, a_10_2_2_2_2, COMP_CPP) );
181  INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<double>(a_10_2_2, a_9_2_2_2_2, a_10_2_2_2_2, COMP_CPP) );
182  INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<double>(a_10_2_2, a_10_2_2_2_2, a_10_2_3_2_2, COMP_CPP) );
183  INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<double>(a_10_2_2, a_10_2_2_2_2, a_10_2_2_3_2, COMP_CPP) );
184  INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<double>(a_10_2_2, a_10_2_2_2_2, a_10_2_2_2_3, COMP_CPP) );
185  INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<double>(a_9_2_2, a_10_2_2_2_2, a_10_2_2_2_2, COMP_CPP) );
186  INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<double>(a_10_3_2, a_10_2_2_2_2, a_10_2_2_2_2, COMP_CPP) );
187  INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<double>(a_10_2_2, a_10_2_2_2_2, a_10_3_2_2_2, COMP_CPP) );
188  INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<double>(a_10_2_3, a_10_2_2_2_2, a_10_3_2_2_2, COMP_ENGINE_MAX) );
189  INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<double>(a_10_2_3, a_10_2_2_2_2, a_10_3_2_2_2, COMP_CPP) );
190 
191 
192  FieldContainer<double> a_9_2(9, 2);
193  FieldContainer<double> a_10_1(10, 1);
194 
195  *outStream << "-> contractDataFieldScalar:\n";
196  INTREPID_TEST_COMMAND( atools.contractDataFieldScalar<double>(a_2_2, a_2_2, a_2_2, COMP_CPP) );
197  INTREPID_TEST_COMMAND( atools.contractDataFieldScalar<double>(a_2_2, a_10_2_2, a_10_2_2, COMP_CPP) );
198  INTREPID_TEST_COMMAND( atools.contractDataFieldScalar<double>(a_10_2_2, a_2_2, a_10_2_2, COMP_CPP) );
199  INTREPID_TEST_COMMAND( atools.contractDataFieldScalar<double>(a_2_2, a_10_2, a_9_2_2, COMP_CPP) );
200  INTREPID_TEST_COMMAND( atools.contractDataFieldScalar<double>(a_2_2, a_10_3, a_10_2_2, COMP_CPP) );
201  INTREPID_TEST_COMMAND( atools.contractDataFieldScalar<double>(a_9_2, a_10_2, a_10_2_2, COMP_CPP) );
202  INTREPID_TEST_COMMAND( atools.contractDataFieldScalar<double>(a_10_2, a_10_2, a_10_3_2, COMP_CPP) );
203  INTREPID_TEST_COMMAND( atools.contractDataFieldScalar<double>(a_10_2, a_10_2, a_10_2_2, COMP_ENGINE_MAX) );
204  INTREPID_TEST_COMMAND( atools.contractDataFieldScalar<double>(a_10_2, a_10_2, a_10_2_2, COMP_CPP) );
205  INTREPID_TEST_COMMAND( atools.contractDataFieldScalar<double>(a_10_2, a_10_1, a_10_2_2, COMP_CPP) );
206 
207 
208  FieldContainer<double> a_10_1_2(10, 1, 2);
209  FieldContainer<double> a_10_1_3(10, 1, 3);
210 
211  *outStream << "-> contractDataFieldVector:\n";
212  INTREPID_TEST_COMMAND( atools.contractDataFieldVector<double>(a_2_2, a_2_2, a_2_2, COMP_CPP) );
213  INTREPID_TEST_COMMAND( atools.contractDataFieldVector<double>(a_2_2, a_2_2, a_10_2_2_2, COMP_CPP) );
214  INTREPID_TEST_COMMAND( atools.contractDataFieldVector<double>(a_10_2_2, a_10_2_2, a_10_2_2_2, COMP_CPP) );
215  INTREPID_TEST_COMMAND( atools.contractDataFieldVector<double>(a_2_2, a_10_2_2, a_9_2_2_2, COMP_CPP) );
216  INTREPID_TEST_COMMAND( atools.contractDataFieldVector<double>(a_2_2, a_10_2_2, a_10_2_3_2, COMP_CPP) );
217  INTREPID_TEST_COMMAND( atools.contractDataFieldVector<double>(a_2_2, a_10_2_2, a_10_2_2_3, COMP_CPP) );
218  INTREPID_TEST_COMMAND( atools.contractDataFieldVector<double>(a_9_2, a_10_2_2, a_10_2_2_2, COMP_CPP) );
219  INTREPID_TEST_COMMAND( atools.contractDataFieldVector<double>(a_10_2, a_10_2_2, a_10_3_2_2, COMP_CPP) );
220  INTREPID_TEST_COMMAND( atools.contractDataFieldVector<double>(a_10_2, a_10_2_2, a_10_2_2_2, COMP_ENGINE_MAX) );
221  INTREPID_TEST_COMMAND( atools.contractDataFieldVector<double>(a_10_2, a_10_2_2, a_10_2_2_2, COMP_CPP) );
222  INTREPID_TEST_COMMAND( atools.contractDataFieldVector<double>(a_10_2, a_10_1_2, a_10_2_2_2, COMP_CPP) );
223 
224 
225  FieldContainer<double> a_10_1_2_2(10, 1, 2, 2);
226 
227  *outStream << "-> contractDataFieldTensor:\n";
228  INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<double>(a_2_2, a_2_2, a_2_2, COMP_CPP) );
229  INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<double>(a_2_2, a_2_2, a_10_2_2_2_2, COMP_CPP) );
230  INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<double>(a_10_2_2, a_10_2_2_2, a_10_2_2_2_2, COMP_CPP) );
231  INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<double>(a_2_2, a_10_2_2_2, a_9_2_2_2_2, COMP_CPP) );
232  INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<double>(a_2_2, a_10_2_2_2, a_10_2_3_2_2, COMP_CPP) );
233  INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<double>(a_2_2, a_10_2_2_2, a_10_2_2_3_2, COMP_CPP) );
234  INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<double>(a_2_2, a_10_2_2_2, a_10_2_2_2_3, COMP_CPP) );
235  INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<double>(a_9_2, a_10_2_2_2, a_10_2_2_2_2, COMP_CPP) );
236  INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<double>(a_10_2, a_10_2_2_2, a_10_3_2_2_2, COMP_CPP) );
237  INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<double>(a_10_2, a_10_2_2_2, a_10_2_2_2_2, COMP_ENGINE_MAX) );
238  INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<double>(a_10_2, a_10_2_2_2, a_10_2_2_2_2, COMP_CPP) );
239  INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<double>(a_10_2, a_10_1_2_2, a_10_2_2_2_2, COMP_CPP) );
240 
241 
242  FieldContainer<double> a_2(2);
243  FieldContainer<double> a_10(10);
244 
245  *outStream << "-> contractDataDataScalar:\n";
246  INTREPID_TEST_COMMAND( atools.contractDataDataScalar<double>(a_2_2, a_10_2_2, a_10_2_2, COMP_CPP) );
247  INTREPID_TEST_COMMAND( atools.contractDataDataScalar<double>(a_2_2, a_10_2, a_10_2_2, COMP_CPP) );
248  INTREPID_TEST_COMMAND( atools.contractDataDataScalar<double>(a_2_2, a_10_2, a_10_2, COMP_CPP) );
249  INTREPID_TEST_COMMAND( atools.contractDataDataScalar<double>(a_2, a_9_2, a_10_2, COMP_CPP) );
250  INTREPID_TEST_COMMAND( atools.contractDataDataScalar<double>(a_2, a_10_2, a_10_3, COMP_CPP) );
251  INTREPID_TEST_COMMAND( atools.contractDataDataScalar<double>(a_2, a_10_2, a_10_2, COMP_CPP) );
252  INTREPID_TEST_COMMAND( atools.contractDataDataScalar<double>(a_10, a_10_2, a_10_2, COMP_ENGINE_MAX) );
253  INTREPID_TEST_COMMAND( atools.contractDataDataScalar<double>(a_10, a_10_2, a_10_2, COMP_CPP) );
254 
255 
256  *outStream << "-> contractDataDataVector:\n";
257  INTREPID_TEST_COMMAND( atools.contractDataDataVector<double>(a_2_2, a_2_2, a_2_2, COMP_CPP) );
258  INTREPID_TEST_COMMAND( atools.contractDataDataVector<double>(a_2_2, a_10_2_2, a_2_2, COMP_CPP) );
259  INTREPID_TEST_COMMAND( atools.contractDataDataVector<double>(a_2_2, a_10_2_2, a_10_2_2, COMP_CPP) );
260  INTREPID_TEST_COMMAND( atools.contractDataDataVector<double>(a_10, a_9_2_2, a_10_2_2, COMP_CPP) );
261  INTREPID_TEST_COMMAND( atools.contractDataDataVector<double>(a_10, a_10_3_2, a_10_2_2, COMP_CPP) );
262  INTREPID_TEST_COMMAND( atools.contractDataDataVector<double>(a_10, a_10_2_3, a_10_2_2, COMP_CPP) );
263  INTREPID_TEST_COMMAND( atools.contractDataDataVector<double>(a_2, a_10_2_2, a_10_2_2, COMP_CPP) );
264  INTREPID_TEST_COMMAND( atools.contractDataDataVector<double>(a_10, a_10_2_2, a_10_2_2, COMP_ENGINE_MAX) );
265  INTREPID_TEST_COMMAND( atools.contractDataDataVector<double>(a_10, a_10_2_2, a_10_2_2, COMP_CPP) );
266 
267 
268  *outStream << "-> contractDataDataTensor:\n";
269  INTREPID_TEST_COMMAND( atools.contractDataDataTensor<double>(a_2_2, a_2_2, a_2_2, COMP_CPP) );
270  INTREPID_TEST_COMMAND( atools.contractDataDataTensor<double>(a_2_2, a_10_2_2_2, a_2_2, COMP_CPP) );
271  INTREPID_TEST_COMMAND( atools.contractDataDataTensor<double>(a_2_2, a_10_2_2_2, a_10_2_2_2, COMP_CPP) );
272  INTREPID_TEST_COMMAND( atools.contractDataDataTensor<double>(a_10, a_9_2_2_2, a_10_2_2_2, COMP_CPP) );
273  INTREPID_TEST_COMMAND( atools.contractDataDataTensor<double>(a_10, a_10_2_2_2, a_10_3_2_2, COMP_CPP) );
274  INTREPID_TEST_COMMAND( atools.contractDataDataTensor<double>(a_10, a_10_2_2_2, a_10_2_3_2, COMP_CPP) );
275  INTREPID_TEST_COMMAND( atools.contractDataDataTensor<double>(a_10, a_10_2_2_2, a_10_2_2_3, COMP_CPP) );
276  INTREPID_TEST_COMMAND( atools.contractDataDataTensor<double>(a_2, a_10_2_2_2, a_10_2_2_2, COMP_CPP) );
277  INTREPID_TEST_COMMAND( atools.contractDataDataTensor<double>(a_10, a_10_2_2_2, a_10_2_2_2, COMP_ENGINE_MAX) );
278  INTREPID_TEST_COMMAND( atools.contractDataDataTensor<double>(a_10, a_10_2_2_2, a_10_2_2_2, COMP_CPP) );
279 
280 
281 #endif
282 
283  }
284  catch (std::logic_error err) {
285  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
286  *outStream << err.what() << '\n';
287  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
288  errorFlag = -1000;
289  };
290 
291 #ifdef HAVE_INTREPID_DEBUG
292  if (Teuchos::TestForException_getThrowNumber() != endThrowNumber)
293  errorFlag++;
294 #endif
295 
296 
297  *outStream \
298  << "\n"
299  << "===============================================================================\n"\
300  << "| TEST 2: correctness of math operations |\n"\
301  << "===============================================================================\n";
302 
303  outStream->precision(20);
304 
305  try {
306  { // start scope
307  *outStream << "\n************ Checking contractFieldFieldScalar ************\n";
308 
309  int c=5, p=9, l=3, r=7;
310 
311  Kokkos::View<double***> in_c_l_p("in_c_l_p",c, l, p);
312  Kokkos::View<double***> in_c_r_p("in_c_r_p",c, r, p);
313  Kokkos::View<double***> out1_c_l_r("out1_c_l_r",c, l, r);
314  Kokkos::View<double***> out2_c_l_r("out2_c_l_r",c, l, r);
315  double zero = INTREPID_TOL*10000.0;
316 
317  // fill with random numbers
318  for (unsigned int i=0; i<in_c_l_p.dimension(0); i++)
319  for (unsigned int j=0; j<in_c_l_p.dimension(1); j++)
320  for (unsigned int k=0; k<in_c_l_p.dimension(2); k++)
321  in_c_l_p(i,j,k) = Teuchos::ScalarTraits<double>::random();
322 
323  for (unsigned int i=0; i<in_c_r_p.dimension(0); i++)
324  for (unsigned int j=0; j<in_c_r_p.dimension(1); j++)
325  for (unsigned int k=0; k<in_c_r_p.dimension(2); k++)
326  in_c_r_p(i,j,k) = Teuchos::ScalarTraits<double>::random();
327 
328 
329  art::contractFieldFieldScalar<double>(out1_c_l_r, in_c_l_p, in_c_r_p, COMP_CPP);
330  art::contractFieldFieldScalar<double>(out2_c_l_r, in_c_l_p, in_c_r_p, COMP_BLAS);
331  rst::subtract(out1_c_l_r, out2_c_l_r);
332  if (rst::vectorNorm(out1_c_l_r, NORM_ONE) > zero) {
333  *outStream << "\n\nINCORRECT contractFieldFieldScalar (1): check COMP_CPP vs. COMP_BLAS; "
334  << " diff-1norm = " << rst::vectorNorm(out1_c_l_r, NORM_ONE) << "\n\n";
335  errorFlag++;
336  }
337  // with sumInto:
338  Kokkos::deep_copy(out1_c_l_r, 2.0);
339  Kokkos::deep_copy(out2_c_l_r, 2.0);
340  art::contractFieldFieldScalar<double>(out1_c_l_r, in_c_l_p, in_c_r_p, COMP_CPP, true);
341  art::contractFieldFieldScalar<double>(out2_c_l_r, in_c_l_p, in_c_r_p, COMP_BLAS, true);
342  rst::subtract(out1_c_l_r, out2_c_l_r);
343  if (rst::vectorNorm(out1_c_l_r, NORM_ONE) > zero) {
344  *outStream << "\n\nINCORRECT contractFieldFieldScalar (1): check COMP_CPP vs. COMP_BLAS; "
345  << " diff-1norm = " << rst::vectorNorm(out1_c_l_r, NORM_ONE) << "\n\n";
346  errorFlag++;
347  }
348  } // end scope
349 
350  { // start scope
351  *outStream << "\n************ Checking contractFieldFieldVector ************\n";
352 
353  int c=5, p=9, l=3, r=7, d=13;
354 
355  Kokkos::View<double****> in_c_l_p_d("in_c_l_p_d",c, l, p, d);
356  Kokkos::View<double****> in_c_r_p_d("in_c_r_p_d",c, r, p, d);
357  Kokkos::View<double***> out1_c_l_r("out1_c_l_r",c, l, r);
358  Kokkos::View<double***> out2_c_l_r("out2_c_l_r",c, l, r);
359  double zero = INTREPID_TOL*10000.0;
360 
361  // fill with random numbers
362  for (unsigned int i=0; i<in_c_l_p_d.dimension(0); i++)
363  for (unsigned int j=0; j<in_c_l_p_d.dimension(1); j++)
364  for (unsigned int k=0; k<in_c_l_p_d.dimension(2); k++)
365  for (unsigned int l=0; l<in_c_l_p_d.dimension(3); l++)
366  in_c_l_p_d(i,j,k,l) = Teuchos::ScalarTraits<double>::random();
367 
368  for (unsigned int i=0; i<in_c_r_p_d.dimension(0); i++)
369  for (unsigned int j=0; j<in_c_r_p_d.dimension(1); j++)
370  for (unsigned int k=0; k<in_c_r_p_d.dimension(2); k++)
371  for (unsigned int l=0; l<in_c_r_p_d.dimension(3); l++)
372  in_c_r_p_d(i,j,k,l) = Teuchos::ScalarTraits<double>::random();
373 
374 
375  art::contractFieldFieldVector<double>(out1_c_l_r, in_c_l_p_d, in_c_r_p_d, COMP_CPP);
376  art::contractFieldFieldVector<double>(out2_c_l_r, in_c_l_p_d, in_c_r_p_d, COMP_BLAS);
377 
378  rst::subtract(out1_c_l_r, out2_c_l_r);
379  if (rst::vectorNorm(out1_c_l_r, NORM_ONE) > zero) {
380  *outStream << "\n\nINCORRECT contractFieldFieldVector (1): check COMP_CPP vs. COMP_BLAS; "
381  << " diff-1norm = " << rst::vectorNorm(out1_c_l_r, NORM_ONE) << "\n\n";
382  errorFlag++;
383  }
384 
385  // with sumInto:
386  Kokkos::deep_copy(out1_c_l_r, 2.0);
387  Kokkos::deep_copy(out2_c_l_r, 2.0);
388  art::contractFieldFieldVector<double>(out1_c_l_r, in_c_l_p_d, in_c_r_p_d, COMP_CPP, true);
389  art::contractFieldFieldVector<double>(out2_c_l_r, in_c_l_p_d, in_c_r_p_d, COMP_BLAS, true);
390 
391  rst::subtract(out1_c_l_r, out2_c_l_r);
392  if (rst::vectorNorm(out1_c_l_r, NORM_ONE) > zero) {
393  *outStream << "\n\nINCORRECT contractFieldFieldVector (1): check COMP_CPP vs. COMP_BLAS; "
394  << " diff-1norm = " << rst::vectorNorm(out1_c_l_r, NORM_ONE) << "\n\n";
395  errorFlag++;
396  }
397  } // end scope
398 
399  { // start scope
400  *outStream << "\n************ Checking contractFieldFieldTensor ************\n";
401 
402  int c=5, p=9, l=3, r=7, d1=13, d2=5;
403 
404  Kokkos::View<double*****> in_c_l_p_d_d("in_c_l_p_d_d", c, l, p, d1, d2);
405  Kokkos::View<double*****> in_c_r_p_d_d("in_c_r_p_d_d", c, r, p, d1, d2);
406  Kokkos::View<double***> out1_c_l_r("out1_c_l_r", c, l, r);
407  Kokkos::View<double***> out2_c_l_r("out2_c_l_r", c, l, r);
408  double zero = INTREPID_TOL*10000.0;
409 
410  // fill with random numbers
411 
412  for (unsigned int i=0; i<in_c_l_p_d_d.dimension(0); i++)
413  for (unsigned int j=0; j<in_c_l_p_d_d.dimension(1); j++)
414  for (unsigned int k=0; k<in_c_l_p_d_d.dimension(2); k++)
415  for (unsigned int l=0; l<in_c_l_p_d_d.dimension(3); l++)
416  for (unsigned int m=0; m<in_c_l_p_d_d.dimension(4); m++)
417  in_c_l_p_d_d(i,j,k,l,m) = Teuchos::ScalarTraits<double>::random();
418 
419  for (unsigned int i=0; i<in_c_r_p_d_d.dimension(0); i++)
420  for (unsigned int j=0; j<in_c_r_p_d_d.dimension(1); j++)
421  for (unsigned int k=0; k<in_c_r_p_d_d.dimension(2); k++)
422  for (unsigned int l=0; l<in_c_r_p_d_d.dimension(3); l++)
423  for (unsigned int m=0; m<in_c_r_p_d_d.dimension(4); m++)
424  in_c_r_p_d_d(i,j,k,l,m) = Teuchos::ScalarTraits<double>::random();
425 
426 
427  art::contractFieldFieldTensor<double>(out1_c_l_r, in_c_l_p_d_d, in_c_r_p_d_d, COMP_CPP);
428  art::contractFieldFieldTensor<double>(out2_c_l_r, in_c_l_p_d_d, in_c_r_p_d_d, COMP_BLAS);
429 
430  rst::subtract(out1_c_l_r, out2_c_l_r);
431  if (rst::vectorNorm(out1_c_l_r, NORM_ONE) > zero) {
432  *outStream << "\n\nINCORRECT contractFieldFieldTensor (1): check COMP_CPP vs. COMP_BLAS; "
433  << " diff-1norm = " << rst::vectorNorm(out1_c_l_r, NORM_ONE) << "\n\n";
434  errorFlag++;
435  }
436 
437  // with sumInto:
438  Kokkos::deep_copy(out1_c_l_r, 2.0);
439  Kokkos::deep_copy(out2_c_l_r, 2.0);
440 
441  art::contractFieldFieldTensor<double>(out1_c_l_r, in_c_l_p_d_d, in_c_r_p_d_d, COMP_CPP, true);
442  art::contractFieldFieldTensor<double>(out2_c_l_r, in_c_l_p_d_d, in_c_r_p_d_d, COMP_BLAS, true);
443 
444  rst::subtract(out1_c_l_r, out2_c_l_r);
445  if (rst::vectorNorm(out1_c_l_r, NORM_ONE) > zero) {
446  *outStream << "\n\nINCORRECT contractFieldFieldTensor (1): check COMP_CPP vs. COMP_BLAS; "
447  << " diff-1norm = " << rst::vectorNorm(out1_c_l_r, NORM_ONE) << "\n\n";
448  errorFlag++;
449  }
450  } // end scope
451 
452  { // start scope
453  *outStream << "\n************ Checking contractDataFieldScalar ************\n";
454 
455  int c=5, p=9, l=7;
456 
457  Kokkos::View<double***> in_c_l_p("in_c_l_p", c, l, p);
458  Kokkos::View<double**> data_c_p("data_c_p", c, p);
459  Kokkos::View<double**> data_c_1("data_c_1", c, 1);
460  Kokkos::View<double**> out1_c_l("out1_c_l", c, l);
461  Kokkos::View<double**> out2_c_l("out2_c_l", c, l);
462  double zero = INTREPID_TOL*10000.0;
463 
464  // fill with random numbers
465  for (unsigned int i=0; i<in_c_l_p.dimension(0); i++)
466  for (unsigned int j=0; j<in_c_l_p.dimension(1); j++)
467  for (unsigned int k=0; k<in_c_l_p.dimension(2); k++)
468  in_c_l_p(i,j,k) = Teuchos::ScalarTraits<double>::random();
469 
470  for (unsigned int i=0; i<data_c_p.dimension(0); i++)
471  for (unsigned int j=0; j<data_c_p.dimension(1); j++)
472  data_c_p(i,j) = Teuchos::ScalarTraits<double>::random();
473 
474  for (unsigned int i=0; i<data_c_1.dimension(0); i++)
475  for (unsigned int j=0; j<data_c_1.dimension(1); j++)
476  data_c_1(i,j) = Teuchos::ScalarTraits<double>::random();
477 
478 
479 
480  // nonconstant data
481  art::contractDataFieldScalar<double>(out1_c_l, data_c_p, in_c_l_p, COMP_CPP);
482  art::contractDataFieldScalar<double>(out2_c_l, data_c_p, in_c_l_p, COMP_BLAS);
483  rst::subtract(out1_c_l, out2_c_l);
484  if (rst::vectorNorm(out1_c_l, NORM_ONE) > zero) {
485  *outStream << "\n\nINCORRECT contractDataFieldScalar (1): check COMP_CPP vs. COMP_BLAS; "
486  << " diff-1norm = " << rst::vectorNorm(out1_c_l, NORM_ONE) << "\n\n";
487  errorFlag++;
488  }
489  // constant data
490  art::contractDataFieldScalar<double>(out1_c_l, data_c_1, in_c_l_p, COMP_CPP);
491  art::contractDataFieldScalar<double>(out2_c_l, data_c_1, in_c_l_p, COMP_BLAS);
492  rst::subtract(out1_c_l, out2_c_l);
493  if (rst::vectorNorm(out1_c_l, NORM_ONE) > zero) {
494  *outStream << "\n\nINCORRECT contractDataFieldScalar (2): check COMP_CPP vs. COMP_BLAS; "
495  << " diff-1norm = " << rst::vectorNorm(out1_c_l, NORM_ONE) << "\n\n";
496  errorFlag++;
497  }
498  // nonconstant data with sumInto
499  Kokkos::deep_copy(out1_c_l, 2.0);
500  Kokkos::deep_copy(out2_c_l, 2.0);
501  art::contractDataFieldScalar<double>(out1_c_l, data_c_p, in_c_l_p, COMP_CPP, true);
502  art::contractDataFieldScalar<double>(out2_c_l, data_c_p, in_c_l_p, COMP_BLAS, true);
503  rst::subtract(out1_c_l, out2_c_l);
504  if (rst::vectorNorm(out1_c_l, NORM_ONE) > zero) {
505  *outStream << "\n\nINCORRECT contractDataFieldScalar (1): check COMP_CPP vs. COMP_BLAS; "
506  << " diff-1norm = " << rst::vectorNorm(out1_c_l, NORM_ONE) << "\n\n";
507  errorFlag++;
508  }
509  } // end scope
510 
511  { // start scope
512  *outStream << "\n************ Checking contractDataFieldVector ************\n";
513 
514  int c=5, p=9, l=7, d=3;
515 
516  Kokkos::View<double****> in_c_l_p_d("in_c_l_p_d", c, l, p, d);
517  Kokkos::View<double***> data_c_p_d("data_c_p_d", c, p, d);
518  Kokkos::View<double***> data_c_1_d("data_c_1_d", c, 1, d);
519  Kokkos::View<double**> out1_c_l("out1_c_l", c, l);
520  Kokkos::View<double**> out2_c_l("out2_c_l", c, l);
521  double zero = INTREPID_TOL*10000.0;
522 
523  // fill with random numbers
524  for (unsigned int i=0; i<in_c_l_p_d.dimension(0); i++)
525  for (unsigned int j=0; j<in_c_l_p_d.dimension(1); j++)
526  for (unsigned int k=0; k<in_c_l_p_d.dimension(2); k++)
527  for (unsigned int l=0; l<in_c_l_p_d.dimension(3); l++)
528  in_c_l_p_d(i,j,k,l) = Teuchos::ScalarTraits<double>::random();
529 
530  for (unsigned int i=0; i<data_c_p_d.dimension(0); i++)
531  for (unsigned int j=0; j<data_c_p_d.dimension(1); j++)
532  for (unsigned int k=0; k<data_c_p_d.dimension(2); k++)
533  data_c_p_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
534 
535  for (unsigned int i=0; i<data_c_1_d.dimension(0); i++)
536  for (unsigned int j=0; j<data_c_1_d.dimension(1); j++)
537  for (unsigned int k=0; k<data_c_1_d.dimension(2); k++)
538  data_c_1_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
539 
540 
541  // nonconstant data
542  art::contractDataFieldVector<double>(out1_c_l, data_c_p_d, in_c_l_p_d, COMP_CPP);
543  art::contractDataFieldVector<double>(out2_c_l, data_c_p_d, in_c_l_p_d, COMP_BLAS);
544  rst::subtract(out1_c_l, out2_c_l);
545  if (rst::vectorNorm(out1_c_l, NORM_ONE) > zero) {
546  *outStream << "\n\nINCORRECT contractDataFieldVector (1): check COMP_CPP vs. COMP_BLAS; "
547  << " diff-1norm = " << rst::vectorNorm(out1_c_l, NORM_ONE) << "\n\n";
548  errorFlag++;
549  }
550  // constant data
551  art::contractDataFieldVector<double>(out1_c_l, data_c_1_d, in_c_l_p_d, COMP_CPP);
552  art::contractDataFieldVector<double>(out2_c_l, data_c_1_d, in_c_l_p_d, COMP_BLAS);
553  rst::subtract(out1_c_l, out2_c_l);
554  if (rst::vectorNorm(out1_c_l, NORM_ONE) > zero) {
555  *outStream << "\n\nINCORRECT contractDataFieldVector (2): check COMP_CPP vs. COMP_BLAS; "
556  << " diff-1norm = " << rst::vectorNorm(out1_c_l, NORM_ONE) << "\n\n";
557  errorFlag++;
558  }
559  // nonconstant data with sumInto
560  Kokkos::deep_copy(out1_c_l, 2.0);
561  Kokkos::deep_copy(out2_c_l, 2.0);
562  art::contractDataFieldVector<double>(out1_c_l, data_c_p_d, in_c_l_p_d, COMP_CPP, true);
563  art::contractDataFieldVector<double>(out2_c_l, data_c_p_d, in_c_l_p_d, COMP_BLAS, true);
564  rst::subtract(out1_c_l, out2_c_l);
565  if (rst::vectorNorm(out1_c_l, NORM_ONE) > zero) {
566  *outStream << "\n\nINCORRECT contractDataFieldVector (1): check COMP_CPP vs. COMP_BLAS; "
567  << " diff-1norm = " << rst::vectorNorm(out1_c_l, NORM_ONE) << "\n\n";
568  errorFlag++;
569  }
570  } // end scope
571 
572  { // start scope
573  *outStream << "\n************ Checking contractDataFieldTensor ************\n";
574 
575  int c=5, p=9, l=7, d1=3, d2=13;
576 
577  Kokkos::View<double*****> in_c_l_p_d_d("in_c_l_p_d_d", c, l, p, d1, d2);
578  Kokkos::View<double****> data_c_p_d_d("data_c_p_d_d", c, p, d1, d2);
579  Kokkos::View<double****> data_c_1_d_d("data_c_1_d_d", c, 1, d1, d2);
580  Kokkos::View<double**> out1_c_l("out1_c_l", c, l);
581  Kokkos::View<double**> out2_c_l("out2_c_l", c, l);
582  double zero = INTREPID_TOL*10000.0;
583 
584  // fill with random numbers
585  for (unsigned int i=0; i<in_c_l_p_d_d.dimension(0); i++)
586  for (unsigned int j=0; j<in_c_l_p_d_d.dimension(1); j++)
587  for (unsigned int k=0; k<in_c_l_p_d_d.dimension(2); k++)
588  for (unsigned int l=0; l<in_c_l_p_d_d.dimension(3); l++)
589  for (unsigned int m=0; m<in_c_l_p_d_d.dimension(4); m++)
590  in_c_l_p_d_d(i,j,k,l,m) = Teuchos::ScalarTraits<double>::random();
591 
592  for (unsigned int i=0; i<data_c_p_d_d.dimension(0); i++)
593  for (unsigned int j=0; j<data_c_p_d_d.dimension(1); j++)
594  for (unsigned int k=0; k<data_c_p_d_d.dimension(2); k++)
595  for (unsigned int l=0; l<data_c_p_d_d.dimension(3); l++)
596  data_c_p_d_d(i,j,k,l) = Teuchos::ScalarTraits<double>::random();
597 
598  for (unsigned int i=0; i<data_c_1_d_d.dimension(0); i++)
599  for (unsigned int j=0; j<data_c_1_d_d.dimension(1); j++)
600  for (unsigned int k=0; k<data_c_1_d_d.dimension(2); k++)
601  for (unsigned int l=0; l<data_c_1_d_d.dimension(3); l++)
602  data_c_1_d_d(i,j,k,l) = Teuchos::ScalarTraits<double>::random();
603 
604 
605  // nonconstant data
606  art::contractDataFieldTensor<double>(out1_c_l, data_c_p_d_d, in_c_l_p_d_d, COMP_CPP);
607  art::contractDataFieldTensor<double>(out2_c_l, data_c_p_d_d, in_c_l_p_d_d, COMP_BLAS);
608  rst::subtract(out1_c_l, out2_c_l);
609  if (rst::vectorNorm(out1_c_l, NORM_ONE) > zero) {
610  *outStream << "\n\nINCORRECT contractDataFieldTensor (1): check COMP_CPP vs. COMP_BLAS; "
611  << " diff-1norm = " << rst::vectorNorm(out1_c_l, NORM_ONE) << "\n\n";
612  errorFlag++;
613  }
614  // constant data
615  art::contractDataFieldTensor<double>(out1_c_l, data_c_1_d_d, in_c_l_p_d_d, COMP_CPP);
616  art::contractDataFieldTensor<double>(out2_c_l, data_c_1_d_d, in_c_l_p_d_d, COMP_BLAS);
617  rst::subtract(out1_c_l, out2_c_l);
618  if (rst::vectorNorm(out1_c_l, NORM_ONE) > zero) {
619  *outStream << "\n\nINCORRECT contractDataFieldTensor (1): check COMP_CPP vs. COMP_BLAS; "
620  << " diff-1norm = " << rst::vectorNorm(out1_c_l, NORM_ONE) << "\n\n";
621  errorFlag++;
622  }
623  // nonconstant data with sumInto
624  Kokkos::deep_copy(out1_c_l, 2.0);
625  Kokkos::deep_copy(out2_c_l, 2.0);
626  art::contractDataFieldTensor<double>(out1_c_l, data_c_p_d_d, in_c_l_p_d_d, COMP_CPP, true);
627  art::contractDataFieldTensor<double>(out2_c_l, data_c_p_d_d, in_c_l_p_d_d, COMP_BLAS, true);
628  rst::subtract(out1_c_l, out2_c_l);
629  if (rst::vectorNorm(out1_c_l, NORM_ONE) > zero) {
630  *outStream << "\n\nINCORRECT contractDataFieldTensor (1): check COMP_CPP vs. COMP_BLAS; "
631  << " diff-1norm = " << rst::vectorNorm(out1_c_l, NORM_ONE) << "\n\n";
632  errorFlag++;
633  }
634  } // end scope
635 
636  { // start scope
637  *outStream << "\n************ Checking contractDataDataScalar ************\n";
638 
639  int c=5, p=9;
640 
641  Kokkos::View<double**> inl_c_p("inl_c_p", c, p);
642  Kokkos::View<double**> inr_c_p("inr_c_p", c, p);
643  Kokkos::View<double*> out1_c("out1_c", c);
644  Kokkos::View<double*> out2_c("out2_c", c);
645  double zero = INTREPID_TOL*10000.0;
646 
647  // fill with random numbers
648  for (unsigned int i=0; i<inl_c_p.dimension(0); i++)
649  for (unsigned int j=0; j<inl_c_p.dimension(1); j++)
650  inl_c_p(i,j) = Teuchos::ScalarTraits<double>::random();
651 
652  for (unsigned int i=0; i<inr_c_p.dimension(0); i++)
653  for (unsigned int j=0; j<inr_c_p.dimension(1); j++)
654  inr_c_p(i,j) = Teuchos::ScalarTraits<double>::random();
655 
656 
657  art::contractDataDataScalar<double>(out1_c, inl_c_p, inr_c_p, COMP_CPP);
658  art::contractDataDataScalar<double>(out2_c, inl_c_p, inr_c_p, COMP_BLAS);
659  rst::subtract(out1_c, out2_c);
660  if (rst::vectorNorm(out1_c, NORM_ONE) > zero) {
661  *outStream << "\n\nINCORRECT contractDataDataScalar (1): check COMP_CPP vs. COMP_BLAS; "
662  << " diff-1norm = " << rst::vectorNorm(out1_c, NORM_ONE) << "\n\n";
663  errorFlag++;
664  }
665  // with sumInto:
666  Kokkos::deep_copy(out1_c, 2.0);
667  Kokkos::deep_copy(out2_c, 2.0);
668  art::contractDataDataScalar<double>(out1_c, inl_c_p, inr_c_p, COMP_CPP, true);
669  art::contractDataDataScalar<double>(out2_c, inl_c_p, inr_c_p, COMP_BLAS, true);
670  rst::subtract(out1_c, out2_c);
671  if (rst::vectorNorm(out1_c, NORM_ONE) > zero) {
672  *outStream << "\n\nINCORRECT contractDataDataScalar (1): check COMP_CPP vs. COMP_BLAS; "
673  << " diff-1norm = " << rst::vectorNorm(out1_c, NORM_ONE) << "\n\n";
674  errorFlag++;
675  }
676  } // end scope
677 
678  { // start scope
679  *outStream << "\n************ Checking contractDataDataVector ************\n";
680 
681  int c=5, p=9, d=13;
682 
683  Kokkos::View<double***> inl_c_p_d("inl_c_p_d", c, p, d);
684  Kokkos::View<double***> inr_c_p_d("inr_c_p_d", c, p, d);
685  Kokkos::View<double*> out1_c("out1_c", c);
686  Kokkos::View<double*> out2_c("out2_c", c);
687  double zero = INTREPID_TOL*10000.0;
688 
689  // fill with random numbers
690  for (unsigned int i=0; i<inl_c_p_d.dimension(0); i++)
691  for (unsigned int j=0; j<inl_c_p_d.dimension(1); j++)
692  for (unsigned int k=0; k<inl_c_p_d.dimension(2); k++)
693  inl_c_p_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
694 
695  for (unsigned int i=0; i<inr_c_p_d.dimension(0); i++)
696  for (unsigned int j=0; j<inr_c_p_d.dimension(1); j++)
697  for (unsigned int k=0; k<inr_c_p_d.dimension(2); k++)
698  inr_c_p_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
699 
700 
701  art::contractDataDataVector<double>(out1_c, inl_c_p_d, inr_c_p_d, COMP_CPP);
702  art::contractDataDataVector<double>(out2_c, inl_c_p_d, inr_c_p_d, COMP_BLAS);
703 
704  rst::subtract(out1_c, out2_c);
705  if (rst::vectorNorm(out1_c, NORM_ONE) > zero) {
706  *outStream << "\n\nINCORRECT contractDataDataVector (1): check COMP_CPP vs. COMP_BLAS; "
707  << " diff-1norm = " << rst::vectorNorm(out1_c, NORM_ONE) << "\n\n";
708  errorFlag++;
709  }
710 
711  // with sumInto:
712  Kokkos::deep_copy(out1_c, 2.0);
713  Kokkos::deep_copy(out2_c, 2.0);
714 
715  art::contractDataDataVector<double>(out1_c, inl_c_p_d, inr_c_p_d, COMP_CPP, true);
716  art::contractDataDataVector<double>(out2_c, inl_c_p_d, inr_c_p_d, COMP_BLAS, true);
717 
718  rst::subtract(out1_c, out2_c);
719  if (rst::vectorNorm(out1_c, NORM_ONE) > zero) {
720  *outStream << "\n\nINCORRECT contractDataDataVector (1): check COMP_CPP vs. COMP_BLAS; "
721  << " diff-1norm = " << rst::vectorNorm(out1_c, NORM_ONE) << "\n\n";
722  errorFlag++;
723  }
724  } // end scope
725 
726  { // start scope
727  *outStream << "\n************ Checking contractDataDataTensor ************\n";
728 
729  int c=5, p=9, d1=13, d2=5;
730 
731  Kokkos::View<double****> inl_c_p_d_d("inl_c_p_d_d", c, p, d1, d2);
732  Kokkos::View<double****> inr_c_p_d_d("inr_c_p_d_d", c, p, d1, d2);
733  Kokkos::View<double*> out1_c("out1_c", c);
734  Kokkos::View<double*> out2_c("out2_c", c);
735  double zero = INTREPID_TOL*10000.0;
736 
737  // fill with random numbers
738  for (unsigned int i=0; i<inl_c_p_d_d.dimension(0); i++)
739  for (unsigned int j=0; j<inl_c_p_d_d.dimension(1); j++)
740  for (unsigned int k=0; k<inl_c_p_d_d.dimension(2); k++)
741  for (unsigned int l=0; l<inl_c_p_d_d.dimension(3); l++)
742  inl_c_p_d_d(i,j,k,l) = Teuchos::ScalarTraits<double>::random();
743 
744  for (unsigned int i=0; i<inr_c_p_d_d.dimension(0); i++)
745  for (unsigned int j=0; j<inr_c_p_d_d.dimension(1); j++)
746  for (unsigned int k=0; k<inr_c_p_d_d.dimension(2); k++)
747  for (unsigned int l=0; l<inr_c_p_d_d.dimension(3); l++)
748  inr_c_p_d_d(i,j,k,l) = Teuchos::ScalarTraits<double>::random();
749 
750  art::contractDataDataTensor<double>(out1_c, inl_c_p_d_d, inr_c_p_d_d, COMP_CPP);
751  art::contractDataDataTensor<double>(out2_c, inl_c_p_d_d, inr_c_p_d_d, COMP_BLAS);
752 
753  rst::subtract(out1_c, out2_c);
754  if (rst::vectorNorm(out1_c, NORM_ONE) > zero) {
755  *outStream << "\n\nINCORRECT contractDataDataTensor (1): check COMP_CPP vs. COMP_BLAS; "
756  << " diff-1norm = " << rst::vectorNorm(out1_c, NORM_ONE) << "\n\n";
757  errorFlag++;
758  }
759 
760  // with sumInto:
761  Kokkos::deep_copy(out1_c, 2.0);
762  Kokkos::deep_copy(out2_c, 2.0);
763 
764  art::contractDataDataTensor<double>(out1_c, inl_c_p_d_d, inr_c_p_d_d, COMP_CPP, true);
765  art::contractDataDataTensor<double>(out2_c, inl_c_p_d_d, inr_c_p_d_d, COMP_BLAS, true);
766 
767  rst::subtract(out1_c, out2_c);
768  if (rst::vectorNorm(out1_c, NORM_ONE) > zero) {
769  *outStream << "\n\nINCORRECT contractDataDataTensor (1): check COMP_CPP vs. COMP_BLAS; "
770  << " diff-1norm = " << rst::vectorNorm(out1_c, NORM_ONE) << "\n\n";
771  errorFlag++;
772  }
773  } // end scope
774 
775  /******************************************/
776  *outStream << "\n";
777  }
778  catch (std::logic_error err) {
779  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
780  *outStream << err.what() << '\n';
781  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
782  errorFlag = -1000;
783  };
784 
785 
786  if (errorFlag != 0)
787  std::cout << "End Result: TEST FAILED\n";
788  else
789  std::cout << "End Result: TEST PASSED\n";
790 
791  // reset format state of std::cout
792  std::cout.copyfmt(oldFormatState);
793  Kokkos::finalize();
794  return errorFlag;
795 }
Intrepid_ArrayTools.hpp
Header file for utility class to provide array tools, such as tensor contractions,...
Intrepid_FieldContainer.hpp
Header file for utility class to provide multidimensional containers.
Intrepid_RealSpaceTools.hpp
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D.
Intrepid::ArrayTools
Utility class that provides methods for higher-order algebraic manipulation of user-defined arrays,...
Definition: Intrepid_ArrayTools.hpp:71
Intrepid::ArrayTools::contractFieldFieldTensor
static void contractFieldFieldTensor(ArrayOutFields &outputFields, const ArrayInFieldsLeft &leftFields, const ArrayInFieldsRight &rightFields, const ECompEngine compEngine, const bool sumInto=false)
Contracts the "point" and "space" dimensions P, D1, and D2 of two rank-5 containers with dimensions (...
Definition: Intrepid_ArrayToolsDefContractions.hpp:189
Intrepid::ArrayTools::contractDataDataScalar
static void contractDataDataScalar(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight, const ECompEngine compEngine, const bool sumInto=false)
Contracts the "point" dimensions P of rank-2 containers with dimensions (C,P), and returns the result...
Definition: Intrepid_ArrayToolsDefContractions.hpp:557
Intrepid::ArrayTools::contractDataFieldScalar
static void contractDataFieldScalar(ArrayOutFields &outputFields, const ArrayInData &inputData, const ArrayInFields &inputFields, const ECompEngine compEngine, const bool sumInto=false)
Contracts the "point" dimensions P of a rank-3 containers and a rank-2 container with dimensions (C,...
Definition: Intrepid_ArrayToolsDefContractions.hpp:267
Intrepid::ArrayTools::contractDataDataVector
static void contractDataDataVector(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight, const ECompEngine compEngine, const bool sumInto=false)
Contracts the "point" and "space" dimensions P and D of rank-3 containers with dimensions (C,...
Definition: Intrepid_ArrayToolsDefContractions.hpp:607
Intrepid::ArrayTools::contractFieldFieldVector
static void contractFieldFieldVector(ArrayOutFields &outputFields, const ArrayInFieldsLeft &leftFields, const ArrayInFieldsRight &rightFields, const ECompEngine compEngine, const bool sumInto=false)
Contracts the "point" and "space" dimensions P and D1 of two rank-4 containers with dimensions (C,...
Definition: Intrepid_ArrayToolsDefContractions.hpp:117
Intrepid::ArrayTools::contractDataFieldVector
static void contractDataFieldVector(ArrayOutFields &outputFields, const ArrayInData &inputData, const ArrayInFields &inputFields, const ECompEngine compEngine, const bool sumInto=false)
Contracts the "point" and "space" dimensions P and D of a rank-4 container and a rank-3 container wit...
Definition: Intrepid_ArrayToolsDefContractions.hpp:354
Intrepid::ArrayTools::contractDataDataTensor
static void contractDataDataTensor(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight, const ECompEngine compEngine, const bool sumInto=false)
Contracts the "point" and "space" dimensions P, D1 and D2 of rank-4 containers with dimensions (C,...
Definition: Intrepid_ArrayToolsDefContractions.hpp:665
Intrepid::ArrayTools::contractFieldFieldScalar
static void contractFieldFieldScalar(ArrayOutFields &outputFields, const ArrayInFieldsLeft &leftFields, const ArrayInFieldsRight &rightFields, const ECompEngine compEngine, const bool sumInto=false)
Contracts the "point" dimension P of two rank-3 containers with dimensions (C,L,P) and (C,...
Definition: Intrepid_ArrayToolsDefContractions.hpp:51
Intrepid::ArrayTools::contractDataFieldTensor
static void contractDataFieldTensor(ArrayOutFields &outputFields, const ArrayInData &inputData, const ArrayInFields &inputFields, const ECompEngine compEngine, const bool sumInto=false)
Contracts the "point" and "space" dimensions P, D1 and D2 of a rank-5 container and a rank-4 containe...
Definition: Intrepid_ArrayToolsDefContractions.hpp:450
Intrepid::RealSpaceTools
Implementation of basic linear algebra functionality in Euclidean space.
Definition: Intrepid_RealSpaceTools.hpp:68
Generated by doxygen 1.9.1