Intrepid
test_03_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: dot multiply |\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 #ifdef HAVE_INTREPID_DEBUG
110  int beginThrowNumber = Teuchos::TestForException_getThrowNumber();
111  int endThrowNumber = beginThrowNumber + 36;
112 #endif
113  typedef ArrayTools art;
114  typedef RealSpaceTools<double> rst;
115 #ifdef HAVE_INTREPID_DEBUG
116  ArrayTools atools;
117 #endif
118  *outStream \
119  << "\n"
120  << "===============================================================================\n"\
121  << "| TEST 1: exceptions |\n"\
122  << "===============================================================================\n";
123 
124  try{
125 
126 #ifdef HAVE_INTREPID_DEBUG
127  FieldContainer<double> a_2(2);
128  FieldContainer<double> a_2_2(2, 2);
129  FieldContainer<double> a_3_2(3, 2);
130  FieldContainer<double> a_2_2_2(2, 2, 2);
131  FieldContainer<double> a_2_2_2_2(2, 2, 2, 2);
132  FieldContainer<double> a_10_1(10, 1);
133  FieldContainer<double> a_10_2(10, 2);
134  FieldContainer<double> a_10_3(10, 3);
135  FieldContainer<double> a_10_1_2(10, 1, 2);
136  FieldContainer<double> a_10_2_2(10, 2, 2);
137  FieldContainer<double> a_10_2_2_2(10, 2, 2, 2);
138  FieldContainer<double> a_10_2_2_2_2(10, 2, 2, 2, 2);
139  FieldContainer<double> a_9_2_2(9, 2, 2);
140  FieldContainer<double> a_10_3_2(10, 3, 2);
141  FieldContainer<double> a_10_2_3(10, 2, 3);
142  FieldContainer<double> a_10_2_2_3(10, 2, 2, 3);
143  FieldContainer<double> a_10_1_2_2(10, 1, 2, 2);
144 
145  *outStream << "-> dotMultiplyDataField:\n";
146  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_2_2, a_2, a_2_2) );
147  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_2_2, a_10_2_2, a_10_2_2_2_2) );
148  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2, a_10_2_2, a_10_2_2_2) );
149  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_9_2_2, a_10_2_2_2) );
150  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_10_3_2, a_10_2_2_2) );
151  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_10_2_3, a_10_2_2_2) );
152  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_10_2_2_3, a_10_2_2_2_2) );
153  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_2_2_2, a_10_2_2_2, a_10_2_2_2_2) );
154  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_3_2, a_10_2_2_2, a_10_2_2_2_2) );
155  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_3, a_10_2_2_2, a_10_2_2_2_2) );
156  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_10_2_2_2, a_10_2_2_2_2) );
157  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_10_1_2_2, a_10_2_2_2_2) );
158  //
159  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_2_2, a_2, a_2) );
160  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_2_2, a_10_2_2, a_10_2) );
161  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2, a_10_2_2, a_10_2_2) );
162  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_10_2_2, a_10_3_2) );
163  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_10_2_2, a_10_2_2) );
164  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_3, a_10_2_2, a_2_2_2) );
165  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_9_2_2, a_2_2_2) );
166  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_10_2_3, a_2_2_2) );
167  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_10_2_2_3, a_2_2_2_2) );
168  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_10_2_2_2, a_2_2_2_2) );
169  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_10_1, a_2_2) );
170 
171  *outStream << "-> dotMultiplyDataData:\n";
172  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_2_2, a_2, a_2_2) );
173  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_2_2, a_10_2_2, a_10_2_2_2) );
174  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2_2, a_10_2_2, a_10_2_2) );
175  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_9_2_2, a_10_2_2) );
176  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_3_2, a_10_2_2) );
177  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_2_3, a_10_2_2) );
178  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_2_2_2, a_10_2_2_3) );
179  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_9_2_2, a_9_2_2) );
180  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_3_2, a_10_3_2) );
181  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_2, a_10_2) );
182  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_2_2, a_10_2_2) );
183  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_2_2_2, a_10_2_2_2) );
184  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_1, a_10_2) );
185  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_1_2, a_10_2_2) );
186  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_1_2_2, a_10_2_2_2) );
187  //
188  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_2_2, a_10_2_2_2_2, a_10_2) );
189  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_2_2, a_10_2_2, a_2) );
190  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2_2, a_10_2_2, a_10_2) );
191  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_2_2, a_10_2) );
192  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_3, a_10_2_2, a_2_2) );
193  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_9_2_2, a_2_2) );
194  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_2_3, a_2_2) );
195  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_2_2_3, a_2_2_2) );
196  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_2, a_2) );
197  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_2_2, a_2_2) );
198  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_2_2_2, a_2_2_2) );
199  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_1, a_2) );
200  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_1_2, a_2_2) );
201  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_1_2_2, a_2_2_2) );
202 #endif
203 
204  }
205  catch (std::logic_error err) {
206  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
207  *outStream << err.what() << '\n';
208  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
209  errorFlag = -1000;
210  };
211 
212 #ifdef HAVE_INTREPID_DEBUG
213  if (Teuchos::TestForException_getThrowNumber() != endThrowNumber)
214  errorFlag++;
215 #endif
216  *outStream \
217  << "\n"
218  << "===============================================================================\n"\
219  << "| TEST 2: correctness of math operations |\n"\
220  << "===============================================================================\n";
221 
222  outStream->precision(20);
223 
224  try {
225  { // start scope
226  *outStream << "\n************ Checking dotMultiplyDataField, (branch 1) ************\n";
227 
228  int c=5, p=9, f=7, d1=6, d2=14;
229 
230  Kokkos::View<double***> in_c_f_p("in_c_f_p", c, f, p);
231  Kokkos::View<double****> in_c_f_p_d("in_c_f_p_d", c, f, p, d1);
232  Kokkos::View<double*****> in_c_f_p_d_d("in_c_f_p_d_d", c, f, p, d1, d2);
233  Kokkos::View<double**> data_c_p("data_c_p", c, p);
234  Kokkos::View<double***> data_c_p_d("data_c_p_d", c, p, d1);
235  Kokkos::View<double****> data_c_p_d_d("data_c_p_d_d", c, p, d1, d2);
236  Kokkos::View<double**> data_c_1("data_c_1", c, 1);
237  Kokkos::View<double***> data_c_1_d("data_c_1_d", c, 1, d1);
238  Kokkos::View<double****> data_c_1_d_d("data_c_1_d_d", c, 1, d1, d2);
239  Kokkos::View<double***> out_c_f_p("out_c_f_p", c, f, p);
240  Kokkos::View<double***> outSM_c_f_p("outSM_c_f_p", c, f, p);
241  Kokkos::View<double***> outDM_c_f_p("outDM_c_f_p", c, f, p);
242  double zero = INTREPID_TOL*10000.0;
243 
244  // fill with random numbers
245  for (unsigned int i=0; i<in_c_f_p.dimension(0); i++)
246  for (unsigned int j=0; j<in_c_f_p.dimension(1); j++)
247  for (unsigned int k=0; k<in_c_f_p.dimension(2); k++)
248  in_c_f_p(i,j,k) = Teuchos::ScalarTraits<double>::random();
249 
250 
251  // fill with alternating 1's and -1's
252  int previous_value=-1;
253  for (unsigned int i=0; i<in_c_f_p_d.dimension(0); i++)
254  for (unsigned int j=0; j<in_c_f_p_d.dimension(1); j++)
255  for (unsigned int k=0; k<in_c_f_p_d.dimension(2); k++)
256  for (unsigned int l=0; l<in_c_f_p_d.dimension(3); l++){
257  if(previous_value==1){
258  in_c_f_p_d(i,j,k,l) = -1;
259  previous_value=-1;
260  }else{
261  in_c_f_p_d(i,j,k,l) = 1;
262  previous_value=1;
263  }
264  }
265  // fill with 1's
266 
267  for (unsigned int i=0; i<in_c_f_p_d_d.dimension(0); i++)
268  for (unsigned int j=0; j<in_c_f_p_d_d.dimension(1); j++)
269  for (unsigned int k=0; k<in_c_f_p_d_d.dimension(2); k++)
270  for (unsigned int l=0; l<in_c_f_p_d_d.dimension(3); l++)
271  for (unsigned int m=0; m<in_c_f_p_d_d.dimension(4); m++)
272  in_c_f_p_d_d(i,j,k,l,m)=1.0;
273 
274  // fill with random numbers
275  for (unsigned int i=0; i<data_c_p.dimension(0); i++)
276  for (unsigned int j=0; j<data_c_p.dimension(1); j++){
277  data_c_p(i,j) = Teuchos::ScalarTraits<double>::random();
278  }
279 
280  // fill with 1's
281  for (unsigned int i=0; i<data_c_1.dimension(0); i++)
282  for (unsigned int j=0; j<data_c_1.dimension(1); j++){
283  data_c_1(i,j) = Teuchos::ScalarTraits<double>::random();
284  }
285 
286  // fill with 1's
287  for (unsigned int i=0; i<data_c_p_d.dimension(0); i++)
288  for (unsigned int j=0; j<data_c_p_d.dimension(1); j++)
289  for (unsigned int k=0; k<data_c_p_d.dimension(2); k++)
290  data_c_p_d(i,j,k) = 1.0;
291 
292 
293  // fill with 1's
294  for (unsigned int i=0; i<data_c_1_d.dimension(0); i++)
295  for (unsigned int j=0; j<data_c_1_d.dimension(1); j++)
296  for (unsigned int k=0; k<data_c_1_d.dimension(2); k++)
297  data_c_1_d(i,j,k) = 1.0;
298 
299  // fill with 1's
300 
301  for (unsigned int i=0; i<data_c_p_d_d.dimension(0); i++)
302  for (unsigned int j=0; j<data_c_p_d_d.dimension(1); j++)
303  for (unsigned int k=0; k<data_c_p_d_d.dimension(2); k++)
304  for (unsigned int l=0; l<data_c_p_d_d.dimension(3); l++)
305  data_c_p_d_d(i,j,k,l)=1.0;
306 
307  // fill with 1's
308 
309  for (unsigned int i=0; i<data_c_1_d_d.dimension(0); i++)
310  for (unsigned int j=0; j<data_c_1_d_d.dimension(1); j++)
311  for (unsigned int k=0; k<data_c_1_d_d.dimension(2); k++)
312  for (unsigned int l=0; l<data_c_1_d_d.dimension(3); l++)
313  data_c_1_d_d(i,j,k,l)=1.0;
314 
315 
316  art::scalarMultiplyDataField<double>(outSM_c_f_p, data_c_p, in_c_f_p);
317  art::dotMultiplyDataField<double>(outDM_c_f_p, data_c_p, in_c_f_p);
318  rst::subtract(out_c_f_p, outSM_c_f_p, outDM_c_f_p);
319  if (rst::vectorNorm(out_c_f_p, NORM_ONE) > zero) {
320  *outStream << "\n\nINCORRECT dotMultiplyDataField (1): check dot multiply for scalars vs. scalar multiply\n\n";
321  errorFlag = -1000;
322  }
323  art::dotMultiplyDataField<double>(out_c_f_p, data_c_p_d, in_c_f_p_d);
324  if (rst::vectorNorm(out_c_f_p, NORM_ONE) > zero) {
325  *outStream << "\n\nINCORRECT dotMultiplyDataField (2): check dot multiply of orthogonal vectors\n\n";
326  errorFlag = -1000;
327  }
328  art::dotMultiplyDataField<double>(out_c_f_p, data_c_p_d_d, in_c_f_p_d_d);
329  if ((rst::vectorNorm(out_c_f_p, NORM_INF) - d1*d2) > zero) {
330  *outStream << "\n\nINCORRECT dotMultiplyDataField (3): check dot multiply for tensors of 1s\n\n";
331  errorFlag = -1000;
332  }
333  art::scalarMultiplyDataField<double>(outSM_c_f_p, data_c_1, in_c_f_p);
334  art::dotMultiplyDataField<double>(outDM_c_f_p, data_c_1, in_c_f_p);
335  rst::subtract(out_c_f_p, outSM_c_f_p, outDM_c_f_p);
336  if (rst::vectorNorm(out_c_f_p, NORM_ONE) > zero) {
337  *outStream << "\n\nINCORRECT dotMultiplyDataField (4): check dot multiply for scalars vs. scalar multiply\n\n";
338  errorFlag = -1000;
339  }
340  art::dotMultiplyDataField<double>(out_c_f_p, data_c_1_d, in_c_f_p_d);
341  if (rst::vectorNorm(out_c_f_p, NORM_ONE) > zero) {
342  *outStream << "\n\nINCORRECT dotMultiplyDataField (5): check dot multiply of orthogonal vectors\n\n";
343  errorFlag = -1000;
344  }
345  art::dotMultiplyDataField<double>(out_c_f_p, data_c_1_d_d, in_c_f_p_d_d);
346  if ((rst::vectorNorm(out_c_f_p, NORM_INF) - d1*d2) > zero) {
347  *outStream << "\n\nINCORRECT dotMultiplyDataField (6): check dot multiply for tensors of 1s\n\n";
348  errorFlag = -1000;
349  }
350  } // end scope
351 
352  { // start scope
353  *outStream << "\n************ Checking dotMultiplyDataField, (branch 2) ************\n";
354 
355  int c=5, p=9, f=7, d1=6, d2=14;
356 
357  Kokkos::View<double**> in_f_p("in_f_p", f, p);
358  Kokkos::View<double***> in_f_p_d("in_f_p_d", f, p, d1);
359  Kokkos::View<double****> in_f_p_d_d("in_f_p_d_d", f, p, d1, d2);
360  Kokkos::View<double**> data_c_p("data_c_p", c, p);
361  Kokkos::View<double***> data_c_p_d("data_c_p_d", c, p, d1);
362  Kokkos::View<double****> data_c_p_d_d("data_c_p_d_d", c, p, d1, d2);
363  Kokkos::View<double**> data_c_1("data_c_1", c, 1);
364  Kokkos::View<double***> data_c_1_d("data_c_1_d", c, 1, d1);
365  Kokkos::View<double****> data_c_1_d_d("data_c_1_d_d", c, 1, d1, d2);
366  Kokkos::View<double***> out_c_f_p("out_c_f_p", c, f, p);
367  Kokkos::View<double***> outSM_c_f_p("outSM_c_f_p", c, f, p);
368  Kokkos::View<double***> outDM_c_f_p("outDM_c_f_p", c, f, p);
369  double zero = INTREPID_TOL*10000.0;
370 
371  // fill with random numbers
372  for (unsigned int i=0; i<in_f_p.dimension(0); i++)
373  for (unsigned int j=0; j<in_f_p.dimension(1); j++){
374  in_f_p(i,j) = Teuchos::ScalarTraits<double>::random();
375  }
376 
377  // fill with alternating 1's and -1's
378  int previous_value=-1;
379  for (unsigned int i=0; i<in_f_p_d.dimension(0); i++)
380  for (unsigned int j=0; j<in_f_p_d.dimension(1); j++)
381  for (unsigned int k=0; k<in_f_p_d.dimension(2); k++){
382  if(previous_value==1){
383  in_f_p_d(i,j,k) = -1;
384  previous_value = -1;
385  }else{
386  in_f_p_d(i,j,k) = 1;
387  previous_value = 1;
388  }
389  }
390 
391  // fill with 1's
392  for (unsigned int i=0; i<in_f_p_d_d.dimension(0); i++)
393  for (unsigned int j=0; j<in_f_p_d_d.dimension(1); j++)
394  for (unsigned int k=0; k<in_f_p_d_d.dimension(2); k++)
395  for (unsigned int l=0; l<in_f_p_d_d.dimension(3); l++)
396  in_f_p_d_d(i,j,k,l) = 1.0;
397 
398  // fill with random numbers
399  for (unsigned int i=0; i<data_c_p.dimension(0); i++)
400  for (unsigned int j=0; j<data_c_p.dimension(1); j++){
401  data_c_p(i,j) = Teuchos::ScalarTraits<double>::random();
402  }
403 
404  // fill with 1's
405  for (unsigned int i=0; i<data_c_1.dimension(0); i++)
406  for (unsigned int j=0; j<data_c_1.dimension(1); j++){
407  data_c_1(i,j) = Teuchos::ScalarTraits<double>::random();
408  }
409 
410  // fill with 1's
411  for (unsigned int i=0; i<data_c_p_d.dimension(0); i++)
412  for (unsigned int j=0; j<data_c_p_d.dimension(1); j++)
413  for (unsigned int k=0; k<data_c_p_d.dimension(2); k++)
414  data_c_p_d(i,j,k) = 1.0;
415 
416  // fill with 1's
417  for (unsigned int i=0; i<data_c_1_d.dimension(0); i++)
418  for (unsigned int j=0; j<data_c_1_d.dimension(1); j++)
419  for (unsigned int k=0; k<data_c_1_d.dimension(2); k++)
420  data_c_1_d(i,j,k) = 1.0;
421 
422  // fill with 1's
423  for (unsigned int i=0; i<data_c_p_d_d.dimension(0); i++)
424  for (unsigned int j=0; j<data_c_p_d_d.dimension(1); j++)
425  for (unsigned int k=0; k<data_c_p_d_d.dimension(2); k++)
426  for (unsigned int l=0; l<data_c_p_d_d.dimension(3); l++)
427  data_c_p_d_d(i,j,k,l)=1.0;
428 
429  // fill with 1's
430  for (unsigned int i=0; i<data_c_1_d_d.dimension(0); i++)
431  for (unsigned int j=0; j<data_c_1_d_d.dimension(1); j++)
432  for (unsigned int k=0; k<data_c_1_d_d.dimension(2); k++)
433  for (unsigned int l=0; l<data_c_1_d_d.dimension(3); l++)
434  data_c_1_d_d(i,j,k,l)=1.0;
435 
436  art::scalarMultiplyDataField<double>(outSM_c_f_p, data_c_p, in_f_p);
437  art::dotMultiplyDataField<double>(outDM_c_f_p, data_c_p, in_f_p);
438  rst::subtract(out_c_f_p, outSM_c_f_p, outDM_c_f_p);
439  if (rst::vectorNorm(out_c_f_p, NORM_ONE) > zero) {
440  *outStream << "\n\nINCORRECT dotMultiplyDataField (7): check dot multiply for scalars vs. scalar multiply\n\n";
441  errorFlag = -1000;
442  }
443  art::dotMultiplyDataField<double>(out_c_f_p, data_c_p_d, in_f_p_d);
444  if (rst::vectorNorm(out_c_f_p, NORM_ONE) > zero) {
445  *outStream << "\n\nINCORRECT dotMultiplyDataField (8): check dot multiply of orthogonal vectors\n\n";
446  errorFlag = -1000;
447  }
448  art::dotMultiplyDataField<double>(out_c_f_p, data_c_p_d_d, in_f_p_d_d);
449  if ((rst::vectorNorm(out_c_f_p, NORM_INF) - d1*d2) > zero) {
450  *outStream << "\n\nINCORRECT dotMultiplyDataField (9): check dot multiply for tensors of 1s\n\n";
451  errorFlag = -1000;
452  }
453  art::scalarMultiplyDataField<double>(outSM_c_f_p, data_c_1, in_f_p);
454  art::dotMultiplyDataField<double>(outDM_c_f_p, data_c_1, in_f_p);
455  rst::subtract(out_c_f_p, outSM_c_f_p, outDM_c_f_p);
456  if (rst::vectorNorm(out_c_f_p, NORM_ONE) > zero) {
457  *outStream << "\n\nINCORRECT dotMultiplyDataField (10): check dot multiply for scalars vs. scalar multiply\n\n";
458  errorFlag = -1000;
459  }
460  art::dotMultiplyDataField<double>(out_c_f_p, data_c_1_d, in_f_p_d);
461  if (rst::vectorNorm(out_c_f_p, NORM_ONE) > zero) {
462  *outStream << "\n\nINCORRECT dotMultiplyDataField (11): check dot multiply of orthogonal vectors\n\n";
463  errorFlag = -1000;
464  }
465  art::dotMultiplyDataField<double>(out_c_f_p, data_c_1_d_d, in_f_p_d_d);
466  if ((rst::vectorNorm(out_c_f_p, NORM_INF) - d1*d2) > zero) {
467  *outStream << "\n\nINCORRECT dotMultiplyDataField (12): check dot multiply for tensors of 1s\n\n";
468  errorFlag = -1000;
469  }
470  } // end scope
471 
472  { // start scope
473  *outStream << "\n************ Checking dotMultiplyDataData, (branch 1) ************\n";
474 
475  int c=5, p=9, d1=6, d2=14;
476 
477  Kokkos::View<double**> in_c_p("in_c_p", c, p);
478  Kokkos::View<double***> in_c_p_d("in_c_p_d", c, p, d1);
479  Kokkos::View<double****> in_c_p_d_d("in_c_p_d_d", c, p, d1, d2);
480  Kokkos::View<double**> data_c_p("data_c_p", c, p);
481  Kokkos::View<double***> data_c_p_d("data_c_p_d", c, p, d1);
482  Kokkos::View<double****> data_c_p_d_d("data_c_p_d_d", c, p, d1, d2);
483  Kokkos::View<double**> data_c_1("data_c_1", c, 1);
484  Kokkos::View<double***> data_c_1_d("data_c_1_d", c, 1, d1);
485  Kokkos::View<double****> data_c_1_d_d("data_c_1_d_d", c, 1, d1, d2);
486  Kokkos::View<double**> out_c_p("out_c_p", c, p);
487  Kokkos::View<double**> outSM_c_p("outSM_c_p", c, p);
488  Kokkos::View<double**> outDM_c_p("outDM_c_p", c, p);
489  double zero = INTREPID_TOL*10000.0;
490 
491  // fill with random numbers
492  for (unsigned int i=0; i<in_c_p.dimension(0); i++)
493  for (unsigned int j=0; j<in_c_p.dimension(1); j++){
494  in_c_p(i,j) = Teuchos::ScalarTraits<double>::random();
495  }
496 
497  // fill with alternating 1's and -1's
498  int previous_value=-1;
499  for (unsigned int i=0; i<in_c_p_d.dimension(0); i++)
500  for (unsigned int j=0; j<in_c_p_d.dimension(1); j++)
501  for (unsigned int k=0; k<in_c_p_d.dimension(2); k++){
502  if(previous_value==1){
503  in_c_p_d(i,j,k) = -1;
504  previous_value = -1;
505  }else{
506  in_c_p_d(i,j,k) = 1;
507  previous_value = 1;
508  }
509  }
510 
511  // fill with 1's
512  for (unsigned int i=0; i<in_c_p_d_d.dimension(0); i++)
513  for (unsigned int j=0; j<in_c_p_d_d.dimension(1); j++)
514  for (unsigned int k=0; k<in_c_p_d_d.dimension(2); k++)
515  for (unsigned int l=0; l<in_c_p_d_d.dimension(3); l++)
516  in_c_p_d_d(i,j,k,l) = 1.0;
517 
518  // fill with random numbers
519  for (unsigned int i=0; i<data_c_p.dimension(0); i++)
520  for (unsigned int j=0; j<data_c_p.dimension(1); j++){
521  data_c_p(i,j) = Teuchos::ScalarTraits<double>::random();
522  }
523  // fill with 1's
524  for (unsigned int i=0; i<data_c_1.dimension(0); i++)
525  for (unsigned int j=0; j<data_c_1.dimension(1); j++){
526  data_c_1(i,j) = Teuchos::ScalarTraits<double>::random();
527  }
528  // fill with 1's
529  for (unsigned int i=0; i<data_c_p_d.dimension(0); i++)
530  for (unsigned int j=0; j<data_c_p_d.dimension(1); j++)
531  for (unsigned int k=0; k<data_c_p_d.dimension(2); k++)
532  data_c_p_d(i,j,k) = 1.0;
533  // fill with 1's
534  for (unsigned int i=0; i<data_c_1_d.dimension(0); i++)
535  for (unsigned int j=0; j<data_c_1_d.dimension(1); j++)
536  for (unsigned int k=0; k<data_c_1_d.dimension(2); k++)
537  data_c_1_d(i,j,k) = 1.0;
538  // fill with 1's
539  for (unsigned int i=0; i<data_c_p_d_d.dimension(0); i++)
540  for (unsigned int j=0; j<data_c_p_d_d.dimension(1); j++)
541  for (unsigned int k=0; k<data_c_p_d_d.dimension(2); k++)
542  for (unsigned int l=0; l<data_c_p_d_d.dimension(3); l++)
543  data_c_p_d_d(i,j,k,l)=1.0;
544 
545  // fill with 1's
546  for (unsigned int i=0; i<data_c_1_d_d.dimension(0); i++)
547  for (unsigned int j=0; j<data_c_1_d_d.dimension(1); j++)
548  for (unsigned int k=0; k<data_c_1_d_d.dimension(2); k++)
549  for (unsigned int l=0; l<data_c_1_d_d.dimension(3); l++)
550  data_c_1_d_d(i,j,k,l)=1.0;
551 
552 
553  art::scalarMultiplyDataData<double>(outSM_c_p, data_c_p, in_c_p);
554  art::dotMultiplyDataData<double>(outDM_c_p, data_c_p, in_c_p);
555  rst::subtract(out_c_p, outSM_c_p, outDM_c_p);
556  if (rst::vectorNorm(out_c_p, NORM_ONE) > zero) {
557  *outStream << "\n\nINCORRECT dotMultiplyDataData (1): check dot multiply for scalars vs. scalar multiply\n\n";
558  errorFlag = -1000;
559  }
560  art::dotMultiplyDataData<double>(out_c_p, data_c_p_d, in_c_p_d);
561  if (rst::vectorNorm(out_c_p, NORM_ONE) > zero) {
562  *outStream << "\n\nINCORRECT dotMultiplyDataData (2): check dot multiply of orthogonal vectors\n\n";
563  errorFlag = -1000;
564  }
565  art::dotMultiplyDataData<double>(out_c_p, data_c_p_d_d, in_c_p_d_d);
566  if ((rst::vectorNorm(out_c_p, NORM_INF) - d1*d2) > zero) {
567  *outStream << "\n\nINCORRECT dotMultiplyDataData (3): check dot multiply for tensors of 1s\n\n";
568  errorFlag = -1000;
569  }
570  art::scalarMultiplyDataData<double>(outSM_c_p, data_c_1, in_c_p);
571  art::dotMultiplyDataData<double>(outDM_c_p, data_c_1, in_c_p);
572  rst::subtract(out_c_p, outSM_c_p, outDM_c_p);
573  if (rst::vectorNorm(out_c_p, NORM_ONE) > zero) {
574  *outStream << "\n\nINCORRECT dotMultiplyDataData (4): check dot multiply for scalars vs. scalar multiply\n\n";
575  errorFlag = -1000;
576  }
577  art::dotMultiplyDataData<double>(out_c_p, data_c_1_d, in_c_p_d);
578  if (rst::vectorNorm(out_c_p, NORM_ONE) > zero) {
579  *outStream << "\n\nINCORRECT dotMultiplyDataData (5): check dot multiply of orthogonal vectors\n\n";
580  errorFlag = -1000;
581  }
582  art::dotMultiplyDataData<double>(out_c_p, data_c_1_d_d, in_c_p_d_d);
583  if ((rst::vectorNorm(out_c_p, NORM_INF) - d1*d2) > zero) {
584  *outStream << "\n\nINCORRECT dotMultiplyDataData (6): check dot multiply for tensors of 1s\n\n";
585  errorFlag = -1000;
586  }
587  } // end scope
588 
589  { // start scope
590  *outStream << "\n************ Checking dotMultiplyDataData, (branch 2) ************\n";
591 
592  int c=5, p=9, d1=6, d2=14;
593 
594  Kokkos::View<double*> in_p("in_p", p);
595  Kokkos::View<double**> in_p_d("in_p_d", p, d1);
596  Kokkos::View<double***> in_p_d_d("in_p_d_d", p, d1, d2);
597  Kokkos::View<double**> data_c_p("data_c_p", c, p);
598  Kokkos::View<double***> data_c_p_d("data_c_p_d", c, p, d1);
599  Kokkos::View<double****> data_c_p_d_d("data_c_p_d_d", c, p, d1, d2);
600  Kokkos::View<double**> data_c_1("data_c_1", c, 1);
601  Kokkos::View<double***> data_c_1_d("data_c_1_d", c, 1, d1);
602  Kokkos::View<double****> data_c_1_d_d("data_c_1_d_d", c, 1, d1, d2);
603  Kokkos::View<double**> out_c_p("out_c_p", c, p);
604  Kokkos::View<double**> outSM_c_p("outSM_c_p", c, p);
605  Kokkos::View<double**> outDM_c_p("outDM_c_p", c, p);
606  double zero = INTREPID_TOL*10000.0;
607 
608  // fill with random numbers
609  for (unsigned int i=0; i<in_p.dimension(0); i++){
610  in_p(i) = Teuchos::ScalarTraits<double>::random();
611  }
612 
613  // fill with alternating 1's and -1's
614  int previous_value=-1;
615  for (unsigned int i=0; i<in_p_d.dimension(0); i++)
616  for (unsigned int j=0; j<in_p_d.dimension(1); j++){
617  if(previous_value==1){
618  in_p_d(i,j) = -1;
619  previous_value = -1;
620  }else{
621  in_p_d(i,j) = 1;
622  previous_value = 1;
623  }
624  }
625 
626  // fill with 1's
627  for (unsigned int i=0; i<in_p_d_d.dimension(0); i++)
628  for (unsigned int j=0; j<in_p_d_d.dimension(1); j++)
629  for (unsigned int k=0; k<in_p_d_d.dimension(2); k++){
630  in_p_d_d(i,j,k) = 1.0;
631  }
632 
633  // fill with random numbers
634  for (unsigned int i=0; i<data_c_p.dimension(0); i++)
635  for (unsigned int j=0; j<data_c_p.dimension(1); j++){
636  data_c_p(i,j) = Teuchos::ScalarTraits<double>::random();
637  }
638  // fill with 1's
639  for (unsigned int i=0; i<data_c_1.dimension(0); i++)
640  for (unsigned int j=0; j<data_c_1.dimension(1); j++){
641  data_c_1(i,j) = Teuchos::ScalarTraits<double>::random();
642  }
643 
644  // fill with 1's
645  for (unsigned int i=0; i<data_c_p_d.dimension(0); i++)
646  for (unsigned int j=0; j<data_c_p_d.dimension(1); j++)
647  for (unsigned int k=0; k<data_c_p_d.dimension(2); k++)
648  data_c_p_d(i,j,k) = 1.0;
649 
650  // fill with 1's
651  for (unsigned int i=0; i<data_c_1_d.dimension(0); i++)
652  for (unsigned int j=0; j<data_c_1_d.dimension(1); j++)
653  for (unsigned int k=0; k<data_c_1_d.dimension(2); k++)
654  data_c_1_d(i,j,k) = 1.0;
655  // fill with 1's
656  for (unsigned int i=0; i<data_c_p_d_d.dimension(0); i++)
657  for (unsigned int j=0; j<data_c_p_d_d.dimension(1); j++)
658  for (unsigned int k=0; k<data_c_p_d_d.dimension(2); k++)
659  for (unsigned int l=0; l<data_c_p_d_d.dimension(3); l++)
660  data_c_p_d_d(i,j,k,l)=1.0;
661  // fill with 1's
662  for (unsigned int i=0; i<data_c_1_d_d.dimension(0); i++)
663  for (unsigned int j=0; j<data_c_1_d_d.dimension(1); j++)
664  for (unsigned int k=0; k<data_c_1_d_d.dimension(2); k++)
665  for (unsigned int l=0; l<data_c_1_d_d.dimension(3); l++)
666  data_c_1_d_d(i,j,k,l)=1.0;
667 
668 
669  art::scalarMultiplyDataData<double>(outSM_c_p, data_c_p, in_p);
670  art::dotMultiplyDataData<double>(outDM_c_p, data_c_p, in_p);
671  rst::subtract(out_c_p, outSM_c_p, outDM_c_p);
672  if (rst::vectorNorm(out_c_p, NORM_ONE) > zero) {
673  *outStream << "\n\nINCORRECT dotMultiplyDataData (7): check dot multiply for scalars vs. scalar multiply\n\n";
674  errorFlag = -1000;
675  }
676  art::dotMultiplyDataData<double>(out_c_p, data_c_p_d, in_p_d);
677  if (rst::vectorNorm(out_c_p, NORM_ONE) > zero) {
678  *outStream << "\n\nINCORRECT dotMultiplyDataData (8): check dot multiply of orthogonal vectors\n\n";
679  errorFlag = -1000;
680  }
681  art::dotMultiplyDataData<double>(out_c_p, data_c_p_d_d, in_p_d_d);
682  if ((rst::vectorNorm(out_c_p, NORM_INF) - d1*d2) > zero) {
683  *outStream << "\n\nINCORRECT dotMultiplyDataData (9): check dot multiply for tensors of 1s\n\n";
684  errorFlag = -1000;
685  }
686  art::scalarMultiplyDataData<double>(outSM_c_p, data_c_1, in_p);
687  art::dotMultiplyDataData<double>(outDM_c_p, data_c_1, in_p);
688  rst::subtract(out_c_p, outSM_c_p, outDM_c_p);
689  if (rst::vectorNorm(out_c_p, NORM_ONE) > zero) {
690  *outStream << "\n\nINCORRECT dotMultiplyDataData (10): check dot multiply for scalars vs. scalar multiply\n\n";
691  errorFlag = -1000;
692  }
693  art::dotMultiplyDataData<double>(out_c_p, data_c_1_d, in_p_d);
694  if (rst::vectorNorm(out_c_p, NORM_ONE) > zero) {
695  *outStream << "\n\nINCORRECT dotMultiplyDataData (11): check dot multiply of orthogonal vectors\n\n";
696  errorFlag = -1000;
697  }
698  art::dotMultiplyDataData<double>(out_c_p, data_c_1_d_d, in_p_d_d);
699  if ((rst::vectorNorm(out_c_p, NORM_INF) - d1*d2) > zero) {
700  *outStream << "\n\nINCORRECT dotMultiplyDataData (12): check dot multiply for tensors of 1s\n\n";
701  errorFlag = -1000;
702  }
703  } // end scope
704  /******************************************/
705  *outStream << "\n";
706  }
707  catch (std::logic_error err) {
708  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
709  *outStream << err.what() << '\n';
710  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
711  errorFlag = -1000;
712  };
713 
714 
715  if (errorFlag != 0)
716  std::cout << "End Result: TEST FAILED\n";
717  else
718  std::cout << "End Result: TEST PASSED\n";
719 
720  // reset format state of std::cout
721  std::cout.copyfmt(oldFormatState);
722  Kokkos::finalize();
723  return errorFlag;
724 }
Intrepid::ArrayTools
Utility class that provides methods for higher-order algebraic manipulation of user-defined arrays,...
Definition: Intrepid_ArrayTools.hpp:71
Intrepid_FieldContainer.hpp
Header file for utility class to provide multidimensional containers.
Intrepid::ArrayTools::dotMultiplyDataField
static void dotMultiplyDataField(ArrayOutFields &outputFields, const ArrayInData &inputDataLeft, const ArrayInFields &inputFields)
There are two use cases: (1) dot product of a rank-3, 4 or 5 container inputFields with dimensions (C...
Definition: Intrepid_ArrayToolsDefDot.hpp:52
Intrepid_RealSpaceTools.hpp
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D.
Intrepid::RealSpaceTools
Implementation of basic linear algebra functionality in Euclidean space.
Definition: Intrepid_RealSpaceTools.hpp:68
Intrepid::FieldContainer< double >
Intrepid::ArrayTools::dotMultiplyDataData
static void dotMultiplyDataData(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight)
There are two use cases: (1) dot product of a rank-2, 3 or 4 container inputDataRight with dimensions...
Definition: Intrepid_ArrayToolsDefDot.hpp:351
Intrepid_ArrayTools.hpp
Header file for utility class to provide array tools, such as tensor contractions,...
Generated by   doxygen 1.8.16