Intrepid
test_02_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"
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 using namespace std;
58 using namespace Intrepid;
59 
60 #define INTREPID_TEST_COMMAND( S ) \
61 { \
62  try { \
63  S ; \
64  } \
65  catch (std::logic_error err) { \
66  *outStream << "Expected Error ----------------------------------------------------------------\n"; \
67  *outStream << err.what() << '\n'; \
68  *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
69  }; \
70 }
71 
72 
73 int main(int argc, char *argv[]) {
74 
75  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
76  Kokkos::initialize();
77  // This little trick lets us print to std::cout only if
78  // a (dummy) command-line argument is provided.
79  int iprint = argc - 1;
80  Teuchos::RCP<std::ostream> outStream;
81  Teuchos::oblackholestream bhs; // outputs nothing
82  if (iprint > 0)
83  outStream = Teuchos::rcp(&std::cout, false);
84  else
85  outStream = Teuchos::rcp(&bhs, false);
86 
87  // Save the format state of the original std::cout.
88  Teuchos::oblackholestream oldFormatState;
89  oldFormatState.copyfmt(std::cout);
90 
91  *outStream \
92  << "===============================================================================\n" \
93  << "| |\n" \
94  << "| Unit Test (ArrayTools) |\n" \
95  << "| |\n" \
96  << "| 1) Array operations: scalar multiply |\n" \
97  << "| |\n" \
98  << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \
99  << "| Denis Ridzal (dridzal@sandia.gov). |\n" \
100  << "| |\n" \
101  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
102  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
103  << "| |\n" \
104  << "===============================================================================\n";
105 
106 
107  int errorFlag = 0;
108 #ifdef HAVE_INTREPID_DEBUG
109  int beginThrowNumber = Teuchos::TestForException_getThrowNumber();
110  int endThrowNumber = beginThrowNumber + 36;
111 #endif
112  typedef ArrayTools art;
113  typedef RealSpaceTools<double> rst;
114 #ifdef HAVE_INTREPID_DEBUG
115  ArrayTools atools;
116 #endif
117  *outStream \
118  << "\n"
119  << "===============================================================================\n"\
120  << "| TEST 1: exceptions |\n"\
121  << "===============================================================================\n";
122 
123  try{
124 
125 #ifdef HAVE_INTREPID_DEBUG
126  FieldContainer<double> a_2_2(2, 2);
127  FieldContainer<double> a_10_2(10, 2);
128  FieldContainer<double> a_10_3(10, 3);
129  FieldContainer<double> a_10_2_2(10, 2, 2);
130  FieldContainer<double> a_10_2_3(10, 2, 3);
131  FieldContainer<double> a_10_3_2(10, 3, 2);
132  FieldContainer<double> a_9_2_2(9, 2, 2);
133 
134  FieldContainer<double> a_10_2_2_2(10, 2, 2, 2);
135  FieldContainer<double> a_9_2_2_2(9, 2, 2, 2);
136  FieldContainer<double> a_10_3_2_2(10, 3, 2, 2);
137  FieldContainer<double> a_10_2_3_2(10, 2, 3, 2);
138  FieldContainer<double> a_10_2_2_3(10, 2, 2, 3);
139 
140  FieldContainer<double> a_10_2_2_2_2(10, 2, 2, 2, 2);
141  FieldContainer<double> a_9_2_2_2_2(9, 2, 2, 2, 2);
142  FieldContainer<double> a_10_3_2_2_2(10, 3, 2, 2, 2);
143  FieldContainer<double> a_10_2_3_2_2(10, 2, 3, 2, 2);
144  FieldContainer<double> a_10_2_2_3_2(10, 2, 2, 3, 2);
145  FieldContainer<double> a_10_2_2_2_3(10, 2, 2, 2, 3);
146 
147  FieldContainer<double> a_9_2(9, 2);
148  FieldContainer<double> a_10_1(10, 1);
149 
150  FieldContainer<double> a_10_1_2(10, 1, 2);
151  FieldContainer<double> a_10_1_3(10, 1, 3);
152 
153  FieldContainer<double> a_10_1_2_2(10, 1, 2, 2);
154 
155  FieldContainer<double> a_2_3_2_2(2, 3, 2, 2);
156  FieldContainer<double> a_2_2_2_2(2, 2, 2, 2);
157  FieldContainer<double> a_2_10(2, 10);
158  FieldContainer<double> a_2(2);
159 
160  *outStream << "-> scalarMultiplyDataField:\n";
161  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_2_2, a_10_2_2, a_2_2) );
162  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_2_2, a_2_2, a_2_2) );
163  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_2_2, a_2_2, a_10_2_2) );
164  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_2_2, a_2_2, a_10_2_2) );
165  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_2_2, a_10_3, a_10_2_2) );
166  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_9_2_2_2_2, a_10_2, a_10_2_2_2_2) );
167  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_3_2_2_2, a_10_2, a_10_2_2_2_2) );
168  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_2_3_2_2, a_10_2, a_10_2_2_2_2) );
169  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_2_2_3_2, a_10_2, a_10_2_2_2_2) );
170  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_2_2_2_3, a_10_2, a_10_2_2_2_2) );
171  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_2_2_2_2, a_10_2, a_10_2_2_2_2) );
172  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_2_2_2_2, a_10_1, a_10_2_2_2_2) );
173  //
174  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_2_2, a_10_2_2, a_2) );
175  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_2_2_2, a_2_2, a_2) );
176  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_2_2_2, a_2_2, a_10_2) );
177  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_2_2, a_10_2, a_2_10) );
178  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_2_2_2_2, a_9_2, a_2_2_2_2) );
179  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_3_2_2_2, a_10_2, a_2_2_2_2) );
180  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_2_3_2_2, a_10_2, a_2_2_2_2) );
181  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_2_2_3_2, a_10_2, a_2_2_2_2) );
182  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_2_2_2_3, a_10_2, a_2_2_2_2) );
183  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_2_2_2_2, a_10_2, a_2_2_2_2) );
184  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_2_2_2_2, a_10_1, a_2_2_2_2) );
185 
186 
187  FieldContainer<double> a_2_2_2(2, 2, 2);
188 
189  *outStream << "-> scalarMultiplyDataData:\n";
190  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_2_2, a_10_2_2, a_2_2) );
191  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_2, a_2_2, a_2) );
192  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_2_2, a_2_2, a_10_2_2) );
193  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_10_2_2, a_2_2, a_10_2_2) );
194  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_10_2_2, a_10_3, a_10_2_2) );
195  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_9_2_2_2, a_10_2, a_10_2_2_2) );
196  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_10_3_2_2, a_10_2, a_10_2_2_2) );
197  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_10_2_3_2, a_10_2, a_10_2_2_2) );
198  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_10_2_2_3, a_10_2, a_10_2_2_2) );
199  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_10_2_2_2, a_10_2, a_10_2_2_2) );
200  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_10_2_2_2, a_10_1, a_10_2_2_2) );
201  //
202  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_2_2, a_10_2_2, a_2) );
203  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_10_2_2_2_2, a_2_2, a_10_2_2_2) );
204  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_10_2_2_2, a_2_2, a_10_2) );
205  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_10_2_2, a_2_2, a_10_2) );
206  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_10_2_2_2, a_9_2, a_2_2_2) );
207  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_10_3_2_2, a_10_2, a_2_2_2) );
208  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_10_2_3_2, a_10_2, a_2_2_2) );
209  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_10_2_2_3, a_10_2, a_2_2_2) );
210  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_10_2_2_2, a_10_2, a_2_2_2) );
211  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_10_2_2_2, a_10_1, a_2_2_2) );
212 #endif
213 
214  }
215  catch (std::logic_error err) {
216  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
217  *outStream << err.what() << '\n';
218  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
219  errorFlag = -1000;
220  };
221 
222 #ifdef HAVE_INTREPID_DEBUG
223  if (Teuchos::TestForException_getThrowNumber() != endThrowNumber)
224  errorFlag++;
225 #endif
226 
227 
228  *outStream \
229  << "\n"
230  << "===============================================================================\n"\
231  << "| TEST 2: correctness of math operations |\n"\
232  << "===============================================================================\n";
233 
234  outStream->precision(20);
235 
236  try {
237  { // start scope
238  *outStream << "\n************ Checking scalarMultiplyDataField, reciprocal=false (branch 1) ************\n";
239 
240  int c=5, p=9, f=7, d1=7, d2=13;
241 
242  Kokkos::View<double***> in_c_f_p("in_c_f_p", c, f, p);
243  Kokkos::View<double****> in_c_f_p_d("in_c_f_p_d", c, f, p, d1);
244  Kokkos::View<double*****> in_c_f_p_d_d("in_c_f_p_d_d", c, f, p, d1, d2);
245  Kokkos::View<double**> data_c_p("data_c_p", c, p);
246  Kokkos::View<double**> datainv_c_p("datainv_c_p", c, p);
247  Kokkos::View<double**> data_c_1("data_c_1", c, 1);
248  Kokkos::View<double**> datainv_c_1("datainv_c_1", c, 1);
249  Kokkos::View<double***> out_c_f_p("out_c_f_p", c, f, p);
250  Kokkos::View<double***> outi_c_f_p("outi_c_f_p", c, f, p);
251  Kokkos::View<double****> out_c_f_p_d("out_c_f_p_d", c, f, p, d1);
252  Kokkos::View<double****> outi_c_f_p_d("outi_c_f_p_d", c, f, p, d1);
253  Kokkos::View<double*****> out_c_f_p_d_d("out_c_f_p_d_d", c, f, p, d1, d2);
254  Kokkos::View<double*****> outi_c_f_p_d_d("outi_c_f_p_d_d", c, f, p, d1, d2);
255  double zero = INTREPID_TOL*10000.0;
256 
257  // fill with random numbers
258  for (unsigned int i=0; i<in_c_f_p.dimension(0); i++)
259  for (unsigned int j=0; j<in_c_f_p.dimension(1); j++)
260  for (unsigned int k=0; k<in_c_f_p.dimension(2); k++)
261  in_c_f_p(i,j,k) = Teuchos::ScalarTraits<double>::random();
262 
263  for (unsigned int i=0; i<in_c_f_p_d.dimension(0); i++)
264  for (unsigned int j=0; j<in_c_f_p_d.dimension(1); j++)
265  for (unsigned int k=0; k<in_c_f_p_d.dimension(2); k++)
266  for (unsigned int l=0; l<in_c_f_p_d.dimension(3); l++)
267  in_c_f_p_d(i,j,k,l) = Teuchos::ScalarTraits<double>::random();
268 
269  for (unsigned int i=0; i<in_c_f_p_d_d.dimension(0); i++)
270  for (unsigned int j=0; j<in_c_f_p_d_d.dimension(1); j++)
271  for (unsigned int k=0; k<in_c_f_p_d_d.dimension(2); k++)
272  for (unsigned int l=0; l<in_c_f_p_d_d.dimension(3); l++)
273  for (unsigned int m=0; m<in_c_f_p_d_d.dimension(4); m++)
274  in_c_f_p_d_d(i,j,k,l,m) = Teuchos::ScalarTraits<double>::random();
275 
276  for (unsigned int i=0; i<data_c_p.dimension(0); i++)
277  for (unsigned int j=0; j<data_c_p.dimension(1); j++){
278  data_c_p(i,j) = Teuchos::ScalarTraits<double>::random();
279  datainv_c_p(i,j) = 1.0 / data_c_p(i,j);
280  }
281 
282  for (unsigned int i=0; i<data_c_1.dimension(0); i++)
283  for (unsigned int j=0; j<data_c_1.dimension(1); j++){
284  data_c_1(i,j) = Teuchos::ScalarTraits<double>::random();
285  datainv_c_1(i,j) = 1.0 / data_c_1(i,j);
286  }
287 
288 
289 
290  art::scalarMultiplyDataField<double>(out_c_f_p, data_c_p, in_c_f_p);
291  art::scalarMultiplyDataField<double>(outi_c_f_p, datainv_c_p, out_c_f_p);
292  rst::subtract(outi_c_f_p, in_c_f_p);
293  if (rst::vectorNorm(outi_c_f_p, NORM_ONE) > zero) {
294  *outStream << "\n\nINCORRECT scalarMultiplyDataField (1): check scalar inverse property\n\n";
295  errorFlag = -1000;
296  }
297  art::scalarMultiplyDataField<double>(out_c_f_p_d, data_c_p, in_c_f_p_d);
298  art::scalarMultiplyDataField<double>(outi_c_f_p_d, datainv_c_p, out_c_f_p_d);
299  rst::subtract(outi_c_f_p_d, in_c_f_p_d);
300  if (rst::vectorNorm(outi_c_f_p_d, NORM_ONE) > zero) {
301  *outStream << "\n\nINCORRECT scalarMultiplyDataField (2): check scalar inverse property\n\n";
302  errorFlag = -1000;
303  }
304  art::scalarMultiplyDataField<double>(out_c_f_p_d_d, data_c_p, in_c_f_p_d_d);
305  art::scalarMultiplyDataField<double>(outi_c_f_p_d_d, datainv_c_p, out_c_f_p_d_d);
306  rst::subtract(outi_c_f_p_d_d, in_c_f_p_d_d);
307  if (rst::vectorNorm(outi_c_f_p_d_d, NORM_ONE) > zero) {
308  *outStream << "\n\nINCORRECT scalarMultiplyDataField (3): check scalar inverse property\n\n";
309  errorFlag = -1000;
310  }
311  art::scalarMultiplyDataField<double>(out_c_f_p_d_d, data_c_1, in_c_f_p_d_d);
312  art::scalarMultiplyDataField<double>(outi_c_f_p_d_d, datainv_c_1, out_c_f_p_d_d);
313  rst::subtract(outi_c_f_p_d_d, in_c_f_p_d_d);
314  if (rst::vectorNorm(outi_c_f_p_d_d, NORM_ONE) > zero) {
315  *outStream << "\n\nINCORRECT scalarMultiplyDataField (4): check scalar inverse property\n\n";
316  errorFlag = -1000;
317  }
318 
319  // fill with constants
320  for (unsigned int i=0; i<in_c_f_p_d_d.dimension(0); i++)
321  for (unsigned int j=0; j<in_c_f_p_d_d.dimension(1); j++)
322  for (unsigned int k=0; k<in_c_f_p_d_d.dimension(2); k++)
323  for (unsigned int l=0; l<in_c_f_p_d_d.dimension(3); l++)
324  for (unsigned int m=0; m<in_c_f_p_d_d.dimension(4); m++)
325  in_c_f_p_d_d(i,j,k,l,m) = 1.0;
326 
327  for (unsigned int i=0; i<data_c_p.dimension(0); i++)
328  for (unsigned int j=0; j<data_c_p.dimension(1); j++){
329  data_c_p(i,j) = 5.0;
330  }
331 
332  for (unsigned int i=0; i<data_c_1.dimension(0); i++)
333  for (unsigned int j=0; j<data_c_1.dimension(1); j++){
334  data_c_1(i,j) = 5.0;
335  }
336 
337  art::scalarMultiplyDataField<double>(out_c_f_p_d_d, data_c_p, in_c_f_p_d_d);
338  if (std::abs(rst::vectorNorm(out_c_f_p_d_d, NORM_ONE) - data_c_p(0,0)*in_c_f_p_d_d(0,0,0,0,0)*c*p*f*d1*d2) > zero) {
339  *outStream << "\n\nINCORRECT scalarMultiplyDataField (5): check result: "
340  << rst::vectorNorm(out_c_f_p_d_d, NORM_ONE) << " != "
341  << data_c_p(0,0)*in_c_f_p_d_d(0,0,0,0,0)*c*f*p*d1*d2 << "\n\n";
342  errorFlag = -1000;
343  }
344  art::scalarMultiplyDataField<double>(out_c_f_p_d_d, data_c_1, in_c_f_p_d_d);
345  if (std::abs(rst::vectorNorm(out_c_f_p_d_d, NORM_ONE) - data_c_1(0,0)*in_c_f_p_d_d(0,0,0,0,0)*c*p*f*d1*d2) > zero) {
346  *outStream << "\n\nINCORRECT scalarMultiplyDataField (6): check result: "
347  << rst::vectorNorm(out_c_f_p_d_d, NORM_ONE) << " != "
348  << data_c_p(0,0)*in_c_f_p_d_d(0,0,0,0,0)*c*f*p*d1*d2 << "\n\n";
349  errorFlag = -1000;
350  }
351  } // end scope
352 
353  { // start scope
354  *outStream << "\n************ Checking scalarMultiplyDataField, reciprocal=false (branch 2) ************\n";
355 
356  int c=5, p=9, f=7, d1=7, d2=13;
357 
358  Kokkos::View<double**> in_f_p("in_f_p", f, p);
359  Kokkos::View<double***> in_f_p_d("in_f_p_d", f, p, d1);
360  Kokkos::View<double****> in_f_p_d_d("in_f_p_d_d", f, p, d1, d2);
361  Kokkos::View<double***> in_c_f_p("in_c_f_p", c, f, p);
362  Kokkos::View<double****> in_c_f_p_d("in_c_f_p_d", c, f, p, d1);
363  Kokkos::View<double*****> in_c_f_p_d_d("in_c_f_p_d_d", c, f, p, d1, d2);
364  Kokkos::View<double**> data_c_p("data_c_p", c, p);
365  Kokkos::View<double**> datainv_c_p("datainv_c_p", c, p);
366  Kokkos::View<double**> data_c_1("data_c_1", c, 1);
367  Kokkos::View<double**> datainv_c_1("datainv_c_1", c, 1);
368  Kokkos::View<double**> data_c_p_one("data_c_p_one", c, p);
369  Kokkos::View<double**> data_c_1_one("data_c_1_one", c, 1);
370  Kokkos::View<double***> out_c_f_p("out_c_f_p", c, f, p);
371  Kokkos::View<double***> outi_c_f_p("outi_c_f_p", c, f, p);
372  Kokkos::View<double****> out_c_f_p_d("out_c_f_p_d", c, f, p, d1);
373  Kokkos::View<double****> outi_c_f_p_d("outi_c_f_p_d", c, f, p, d1);
374  Kokkos::View<double*****> out_c_f_p_d_d("out_c_f_p_d_d", c, f, p, d1, d2);
375  Kokkos::View<double*****> outi_c_f_p_d_d("outi_c_f_p_d_d", c, f, p, d1, d2);
376  double zero = INTREPID_TOL*10000.0;
377 
378  // fill with random numbers
379  for (unsigned int i=0; i<in_f_p.dimension(0); i++)
380  for (unsigned int j=0; j<in_f_p.dimension(1); j++)
381  in_f_p(i,j) = Teuchos::ScalarTraits<double>::random();
382 
383  for (unsigned int i=0; i<in_f_p_d.dimension(0); i++)
384  for (unsigned int j=0; j<in_f_p_d.dimension(1); j++)
385  for (unsigned int k=0; k<in_f_p_d.dimension(2); k++)
386  in_f_p_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
387 
388  for (unsigned int i=0; i<in_f_p_d_d.dimension(0); i++)
389  for (unsigned int j=0; j<in_f_p_d_d.dimension(1); j++)
390  for (unsigned int k=0; k<in_f_p_d_d.dimension(2); k++)
391  for (unsigned int l=0; l<in_f_p_d_d.dimension(3); l++)
392  in_f_p_d_d(i,j,k,l) = Teuchos::ScalarTraits<double>::random();
393 
394  for (unsigned int i=0; i<data_c_p.dimension(0); i++)
395  for (unsigned int j=0; j<data_c_p.dimension(1); j++){
396  data_c_p(i,j) = Teuchos::ScalarTraits<double>::random();
397  datainv_c_p(i,j) = 1.0 / data_c_p(i,j);
398  data_c_p_one(i,j) = 1.0;
399  }
400 
401  for (unsigned int i=0; i<data_c_1.dimension(0); i++)
402  for (unsigned int j=0; j<data_c_1.dimension(1); j++){
403  data_c_1(i,j) = Teuchos::ScalarTraits<double>::random();
404  datainv_c_1(i,j) = 1.0 / data_c_1(i,j);
405  data_c_1_one(i,j) = 1.0;
406 
407  }
408 
409  art::scalarMultiplyDataField<double>(out_c_f_p, data_c_p, in_f_p);
410  art::scalarMultiplyDataField<double>(outi_c_f_p, datainv_c_p, out_c_f_p);
411  art::scalarMultiplyDataField<double>(in_c_f_p, data_c_p_one, in_f_p);
412  rst::subtract(outi_c_f_p, in_c_f_p);
413  if (rst::vectorNorm(outi_c_f_p, NORM_ONE) > zero) {
414  *outStream << "\n\nINCORRECT scalarMultiplyDataField (1): check scalar inverse property\n\n";
415  errorFlag = -1000;
416  }
417  art::scalarMultiplyDataField<double>(out_c_f_p_d, data_c_p, in_f_p_d);
418  art::scalarMultiplyDataField<double>(outi_c_f_p_d, datainv_c_p, out_c_f_p_d);
419  art::scalarMultiplyDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
420  rst::subtract(outi_c_f_p_d, in_c_f_p_d);
421  if (rst::vectorNorm(outi_c_f_p_d, NORM_ONE) > zero) {
422  *outStream << "\n\nINCORRECT scalarMultiplyDataField (2): check scalar inverse property\n\n";
423  errorFlag = -1000;
424  }
425  art::scalarMultiplyDataField<double>(out_c_f_p_d_d, data_c_p, in_f_p_d_d);
426  art::scalarMultiplyDataField<double>(outi_c_f_p_d_d, datainv_c_p, out_c_f_p_d_d);
427  art::scalarMultiplyDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
428  rst::subtract(outi_c_f_p_d_d, in_c_f_p_d_d);
429  if (rst::vectorNorm(outi_c_f_p_d_d, NORM_ONE) > zero) {
430  *outStream << "\n\nINCORRECT scalarMultiplyDataField (3): check scalar inverse property\n\n";
431  errorFlag = -1000;
432  }
433  art::scalarMultiplyDataField<double>(out_c_f_p_d_d, data_c_1, in_f_p_d_d);
434  art::scalarMultiplyDataField<double>(outi_c_f_p_d_d, datainv_c_1, out_c_f_p_d_d);
435  art::scalarMultiplyDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
436  rst::subtract(outi_c_f_p_d_d, in_c_f_p_d_d);
437  if (rst::vectorNorm(outi_c_f_p_d_d, NORM_ONE) > zero) {
438  *outStream << "\n\nINCORRECT scalarMultiplyDataField (4): check scalar inverse property\n\n";
439  errorFlag = -1000;
440  }
441 
442  // fill with constants
443 
444  for (unsigned int i=0; i<in_f_p_d_d.dimension(0); i++)
445  for (unsigned int j=0; j<in_f_p_d_d.dimension(1); j++)
446  for (unsigned int k=0; k<in_f_p_d_d.dimension(2); k++)
447  for (unsigned int l=0; l<in_f_p_d_d.dimension(3); l++)
448  in_f_p_d_d(i,j,k,l) = 1.0;
449 
450 
451  for (unsigned int i=0; i<data_c_p.dimension(0); i++)
452  for (unsigned int j=0; j<data_c_p.dimension(1); j++)
453  data_c_p(i,j) = 5.0;
454 
455  for (unsigned int i=0; i<data_c_1.dimension(0); i++)
456  for (unsigned int j=0; j<data_c_1.dimension(1); j++)
457  data_c_1(i,j) = 5.0;
458 
459  art::scalarMultiplyDataField<double>(out_c_f_p_d_d, data_c_p, in_f_p_d_d);
460  if (std::abs(rst::vectorNorm(out_c_f_p_d_d, NORM_ONE) - data_c_p(0,0)*in_f_p_d_d(0,0,0,0)*c*p*f*d1*d2) > zero) {
461  *outStream << "\n\nINCORRECT scalarMultiplyDataField (5): check result: "
462  << rst::vectorNorm(out_c_f_p_d_d, NORM_ONE) << " != "
463  << data_c_p(0,0)*in_f_p_d_d(0,0,0,0)*c*f*p*d1*d2 << "\n\n";
464  errorFlag = -1000;
465  }
466  art::scalarMultiplyDataField<double>(out_c_f_p_d_d, data_c_1, in_f_p_d_d);
467  if (std::abs(rst::vectorNorm(out_c_f_p_d_d, NORM_ONE) - data_c_1(0,0)*in_f_p_d_d(0,0,0,0)*c*p*f*d1*d2) > zero) {
468  *outStream << "\n\nINCORRECT scalarMultiplyDataField (6): check result: "
469  << rst::vectorNorm(out_c_f_p_d_d, NORM_ONE) << " != "
470  << data_c_1(0,0)*in_f_p_d_d(0,0,0,0)*c*f*p*d1*d2 << "\n\n";
471  errorFlag = -1000;
472  }
473  } // end scope
474 
475  { // start scope
476  *outStream << "\n************ Checking scalarMultiplyDataField, reciprocal=true, i.e. division (branch 1) ************\n";
477 
478  int c=5, p=9, f=7, d1=7, d2=13;
479 
480  Kokkos::View<double***> in_c_f_p("in_c_f_p", c, f, p);
481  Kokkos::View<double****> in_c_f_p_d("in_c_f_p_d", c, f, p, d1);
482  Kokkos::View<double*****> in_c_f_p_d_d("in_c_f_p_d_d", c, f, p, d1, d2);
483  Kokkos::View<double**> data_c_p("data_c_p", c, p);
484  Kokkos::View<double**> datainv_c_p("datainv_c_p", c, p);
485  Kokkos::View<double**> data_c_1("data_c_1", c, 1);
486  Kokkos::View<double**> datainv_c_1("datainv_c_1", c, 1);
487  Kokkos::View<double***> out_c_f_p("out_c_f_p", c, f, p);
488  Kokkos::View<double***> outi_c_f_p("outi_c_f_p", c, f, p);
489  Kokkos::View<double****> out_c_f_p_d("out_c_f_p_d", c, f, p, d1);
490  Kokkos::View<double****> outi_c_f_p_d("outi_c_f_p_d", c, f, p, d1);
491  Kokkos::View<double*****> out_c_f_p_d_d("out_c_f_p_d_d", c, f, p, d1, d2);
492  Kokkos::View<double*****> outi_c_f_p_d_d("outi_c_f_p_d_d", c, f, p, d1, d2);
493  double zero = INTREPID_TOL*10000.0;
494 
495  // fill with random numbers
496  for (unsigned int i=0; i<in_c_f_p.dimension(0); i++)
497  for (unsigned int j=0; j<in_c_f_p.dimension(1); j++)
498  for (unsigned int k=0; k<in_c_f_p.dimension(2); k++)
499  in_c_f_p(i,j,k) = Teuchos::ScalarTraits<double>::random();
500 
501  for (unsigned int i=0; i<in_c_f_p_d.dimension(0); i++)
502  for (unsigned int j=0; j<in_c_f_p_d.dimension(1); j++)
503  for (unsigned int k=0; k<in_c_f_p_d.dimension(2); k++)
504  for (unsigned int l=0; l<in_c_f_p_d.dimension(3); l++)
505  in_c_f_p_d(i,j,k,l) = Teuchos::ScalarTraits<double>::random();
506 
507  for (unsigned int i=0; i<in_c_f_p_d_d.dimension(0); i++)
508  for (unsigned int j=0; j<in_c_f_p_d_d.dimension(1); j++)
509  for (unsigned int k=0; k<in_c_f_p_d_d.dimension(2); k++)
510  for (unsigned int l=0; l<in_c_f_p_d_d.dimension(3); l++)
511  for (unsigned int m=0; m<in_c_f_p_d_d.dimension(4); m++)
512  in_c_f_p_d_d(i,j,k,l,m) = Teuchos::ScalarTraits<double>::random();
513 
514  for (unsigned int i=0; i<data_c_p.dimension(0); i++)
515  for (unsigned int j=0; j<data_c_p.dimension(1); j++){
516  data_c_p(i,j) = Teuchos::ScalarTraits<double>::random();
517  datainv_c_p(i,j) = 1.0 / data_c_p(i,j);
518 
519 }
520 
521  for (unsigned int i=0; i<data_c_1.dimension(0); i++)
522  for (unsigned int j=0; j<data_c_1.dimension(1); j++){
523  data_c_1(i,j) = Teuchos::ScalarTraits<double>::random();
524  datainv_c_1(i,j) = 1.0 / data_c_1(i,j);
525 
526 }
527 
528 
529  art::scalarMultiplyDataField<double>(out_c_f_p, data_c_p, in_c_f_p, true);
530  art::scalarMultiplyDataField<double>(outi_c_f_p, datainv_c_p, out_c_f_p, true);
531  rst::subtract(outi_c_f_p, in_c_f_p);
532  if (rst::vectorNorm(outi_c_f_p, NORM_ONE) > zero) {
533  *outStream << "\n\nINCORRECT scalarMultiplyDataField (1): check scalar inverse property\n\n";
534  errorFlag = -1000;
535  }
536  art::scalarMultiplyDataField<double>(out_c_f_p_d, data_c_p, in_c_f_p_d, true);
537  art::scalarMultiplyDataField<double>(outi_c_f_p_d, datainv_c_p, out_c_f_p_d, true);
538  rst::subtract(outi_c_f_p_d, in_c_f_p_d);
539  if (rst::vectorNorm(outi_c_f_p_d, NORM_ONE) > zero) {
540  *outStream << "\n\nINCORRECT scalarMultiplyDataField (2): check scalar inverse property\n\n";
541  errorFlag = -1000;
542  }
543  art::scalarMultiplyDataField<double>(out_c_f_p_d_d, data_c_p, in_c_f_p_d_d, true);
544  art::scalarMultiplyDataField<double>(outi_c_f_p_d_d, datainv_c_p, out_c_f_p_d_d, true);
545  rst::subtract(outi_c_f_p_d_d, in_c_f_p_d_d);
546  if (rst::vectorNorm(outi_c_f_p_d_d, NORM_ONE) > zero) {
547  *outStream << "\n\nINCORRECT scalarMultiplyDataField (3): check scalar inverse property\n\n";
548  errorFlag = -1000;
549  }
550  art::scalarMultiplyDataField<double>(out_c_f_p_d_d, data_c_1, in_c_f_p_d_d, true);
551  art::scalarMultiplyDataField<double>(outi_c_f_p_d_d, datainv_c_1, out_c_f_p_d_d, true);
552  rst::subtract(outi_c_f_p_d_d, in_c_f_p_d_d);
553  if (rst::vectorNorm(outi_c_f_p_d_d, NORM_ONE) > zero) {
554  *outStream << "\n\nINCORRECT scalarMultiplyDataField (4): check scalar inverse property\n\n";
555  errorFlag = -1000;
556  }
557 
558  // fill with constants
559  for (unsigned int i=0; i<in_c_f_p_d_d.dimension(0); i++)
560  for (unsigned int j=0; j<in_c_f_p_d_d.dimension(1); j++)
561  for (unsigned int k=0; k<in_c_f_p_d_d.dimension(2); k++)
562  for (unsigned int l=0; l<in_c_f_p_d_d.dimension(3); l++)
563  for (unsigned int m=0; m<in_c_f_p_d_d.dimension(4); m++)
564  in_c_f_p_d_d(i,j,k,l,m) = 1.0;
565 
566  for (unsigned int i=0; i<data_c_p.dimension(0); i++)
567  for (unsigned int j=0; j<data_c_p.dimension(1); j++){
568  data_c_p(i,j) = 5.0;
569 
570 
571 }
572 
573  for (unsigned int i=0; i<data_c_1.dimension(0); i++)
574  for (unsigned int j=0; j<data_c_1.dimension(1); j++){
575  data_c_1(i,j) = 5.0;
576 
577 }
578 
579  art::scalarMultiplyDataField<double>(out_c_f_p_d_d, data_c_p, in_c_f_p_d_d, true);
580  if (std::abs(rst::vectorNorm(out_c_f_p_d_d, NORM_ONE) - (1.0/data_c_p(0,0))*in_c_f_p_d_d(0,0,0,0,0)*c*p*f*d1*d2)/rst::vectorNorm(out_c_f_p_d_d, NORM_ONE) > zero) {
581  *outStream << "\n\nINCORRECT scalarMultiplyDataField (5): check result: "
582  << rst::vectorNorm(out_c_f_p_d_d, NORM_ONE) << " != "
583  << (1.0/data_c_p(0,0))*in_c_f_p_d_d(0,0,0,0,0)*c*p*f*d1*d2 << "\n\n";
584  errorFlag = -1000;
585  }
586  art::scalarMultiplyDataField<double>(out_c_f_p_d_d, data_c_1, in_c_f_p_d_d, true);
587  if (std::abs(rst::vectorNorm(out_c_f_p_d_d, NORM_ONE) - (1.0/data_c_1(0,0))*in_c_f_p_d_d(0,0,0,0,0)*c*p*f*d1*d2)/rst::vectorNorm(out_c_f_p_d_d, NORM_ONE) > zero) {
588  *outStream << "\n\nINCORRECT scalarMultiplyDataField (6): check result: "
589  << rst::vectorNorm(out_c_f_p_d_d, NORM_ONE) << " != "
590  << (1.0/data_c_p(0,0))*in_c_f_p_d_d(0,0,0,0,0)*c*p*f*d1*d2 << "\n\n";
591  errorFlag = -1000;
592  }
593 
594  } // end scope
595  { // start scope
596  *outStream << "\n************ Checking scalarMultiplyDataField, reciprocal=true, i.e. division (branch 2) ************\n";
597 
598  int c=5, p=9, f=7, d1=7, d2=13;
599 
600  Kokkos::View<double**> in_f_p("in_f_p", f, p);
601  Kokkos::View<double***> in_f_p_d("in_f_p_d", f, p, d1);
602  Kokkos::View<double****> in_f_p_d_d("in_f_p_d_d", f, p, d1, d2);
603  Kokkos::View<double***> in_c_f_p("in_c_f_p", c, f, p);
604  Kokkos::View<double****> in_c_f_p_d("in_c_f_p_d", c, f, p, d1);
605  Kokkos::View<double*****> in_c_f_p_d_d("in_c_f_p_d_d", c, f, p, d1, d2);
606  Kokkos::View<double**> data_c_p("data_c_p", c, p);
607  Kokkos::View<double**> datainv_c_p("datainv_c_p", c, p);
608  Kokkos::View<double**> data_c_1("data_c_1", c, 1);
609  Kokkos::View<double**> datainv_c_1("datainv_c_1", c, 1);
610  Kokkos::View<double**> data_c_p_one("data_c_p_one", c, p);
611  Kokkos::View<double**> data_c_1_one("data_c_1_one", c, 1);
612  Kokkos::View<double***> out_c_f_p("out_c_f_p", c, f, p);
613  Kokkos::View<double***> outi_c_f_p("outi_c_f_p", c, f, p);
614  Kokkos::View<double****> out_c_f_p_d("out_c_f_p_d", c, f, p, d1);
615  Kokkos::View<double****> outi_c_f_p_d("outi_c_f_p_d", c, f, p, d1);
616  Kokkos::View<double*****> out_c_f_p_d_d("out_c_f_p_d_d", c, f, p, d1, d2);
617  Kokkos::View<double*****> outi_c_f_p_d_d("outi_c_f_p_d_d", c, f, p, d1, d2);
618  double zero = INTREPID_TOL*10000.0;
619 
620  // fill with random numbers
621  for (unsigned int i=0; i<in_f_p.dimension(0); i++)
622  for (unsigned int j=0; j<in_f_p.dimension(1); j++)
623  in_f_p(i,j) = Teuchos::ScalarTraits<double>::random();
624 
625 
626  for (unsigned int i=0; i<in_f_p_d.dimension(0); i++)
627  for (unsigned int j=0; j<in_f_p_d.dimension(1); j++)
628  for (unsigned int k=0; k<in_f_p_d.dimension(2); k++)
629  in_f_p_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
630 
631  for (unsigned int i=0; i<in_f_p_d_d.dimension(0); i++)
632  for (unsigned int j=0; j<in_f_p_d_d.dimension(1); j++)
633  for (unsigned int k=0; k<in_f_p_d_d.dimension(2); k++)
634  for (unsigned int l=0; l<in_f_p_d_d.dimension(3); l++)
635  in_f_p_d_d(i,j,k,l) = Teuchos::ScalarTraits<double>::random();
636 
637  for (unsigned int i=0; i<data_c_p.dimension(0); i++)
638  for (unsigned int j=0; j<data_c_p.dimension(1); j++){
639  data_c_p(i,j) = Teuchos::ScalarTraits<double>::random();
640  datainv_c_p(i,j) = 1.0 / data_c_p(i,j);
641  data_c_p_one(i,j) = 1.0;
642  }
643 
644  for (unsigned int i=0; i<data_c_1.dimension(0); i++)
645  for (unsigned int j=0; j<data_c_1.dimension(1); j++){
646  data_c_1(i,j) = Teuchos::ScalarTraits<double>::random();
647  datainv_c_1(i,j) = 1.0 / data_c_1(i,j);
648  data_c_1_one(i,j) = 1.0;
649  }
650 
651  art::scalarMultiplyDataField<double>(out_c_f_p, data_c_p, in_f_p, true);
652  art::scalarMultiplyDataField<double>(outi_c_f_p, datainv_c_p, out_c_f_p, true);
653  art::scalarMultiplyDataField<double>(in_c_f_p, data_c_p_one, in_f_p);
654  rst::subtract(outi_c_f_p, in_c_f_p);
655  if (rst::vectorNorm(outi_c_f_p, NORM_ONE) > zero) {
656  *outStream << "\n\nINCORRECT scalarMultiplyDataField (1): check scalar inverse property\n\n";
657  errorFlag = -1000;
658  }
659  art::scalarMultiplyDataField<double>(out_c_f_p_d, data_c_p, in_f_p_d, true);
660  art::scalarMultiplyDataField<double>(outi_c_f_p_d, datainv_c_p, out_c_f_p_d, true);
661  art::scalarMultiplyDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
662  rst::subtract(outi_c_f_p_d, in_c_f_p_d);
663  if (rst::vectorNorm(outi_c_f_p_d, NORM_ONE) > zero) {
664  *outStream << "\n\nINCORRECT scalarMultiplyDataField (2): check scalar inverse property\n\n";
665  errorFlag = -1000;
666  }
667  art::scalarMultiplyDataField<double>(out_c_f_p_d_d, data_c_p, in_f_p_d_d, true);
668  art::scalarMultiplyDataField<double>(outi_c_f_p_d_d, datainv_c_p, out_c_f_p_d_d, true);
669  art::scalarMultiplyDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
670  rst::subtract(outi_c_f_p_d_d, in_c_f_p_d_d);
671  if (rst::vectorNorm(outi_c_f_p_d_d, NORM_ONE) > zero) {
672  *outStream << "\n\nINCORRECT scalarMultiplyDataField (3): check scalar inverse property\n\n";
673  errorFlag = -1000;
674  }
675  art::scalarMultiplyDataField<double>(out_c_f_p_d_d, data_c_1, in_f_p_d_d, true);
676  art::scalarMultiplyDataField<double>(outi_c_f_p_d_d, datainv_c_1, out_c_f_p_d_d, true);
677  art::scalarMultiplyDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
678  rst::subtract(outi_c_f_p_d_d, in_c_f_p_d_d);
679  if (rst::vectorNorm(outi_c_f_p_d_d, NORM_ONE) > zero) {
680  *outStream << "\n\nINCORRECT scalarMultiplyDataField (4): check scalar inverse property\n\n";
681  errorFlag = -1000;
682  }
683 
684  // fill with constants
685 
686  for (unsigned int i=0; i<in_f_p_d_d.dimension(0); i++)
687  for (unsigned int j=0; j<in_f_p_d_d.dimension(1); j++)
688  for (unsigned int k=0; k<in_f_p_d_d.dimension(2); k++)
689  for (unsigned int l=0; l<in_f_p_d_d.dimension(3); l++)
690  in_f_p_d_d(i,j,k,l) = 1.0;
691 
692  for (unsigned int i=0; i<data_c_p.dimension(0); i++)
693  for (unsigned int j=0; j<data_c_p.dimension(1); j++){
694  data_c_p(i,j) = 5.0;
695  }
696 
697  for (unsigned int i=0; i<data_c_1.dimension(0); i++)
698  for (unsigned int j=0; j<data_c_1.dimension(1); j++){
699  data_c_1(i,j) = 5.0;
700  }
701 
702  art::scalarMultiplyDataField<double>(out_c_f_p_d_d, data_c_p, in_f_p_d_d, true);
703  if (std::abs(rst::vectorNorm(out_c_f_p_d_d, NORM_ONE) - (1.0/data_c_p(0,0))*in_f_p_d_d(0,0,0,0)*c*p*f*d1*d2)/rst::vectorNorm(out_c_f_p_d_d, NORM_ONE) > zero) {
704  *outStream << "\n\nINCORRECT scalarMultiplyDataField (5): check result: "
705  << rst::vectorNorm(out_c_f_p_d_d, NORM_ONE) << " != "
706  << (1.0/data_c_p(0,0))*in_f_p_d_d(0,0,0,0)*c*p*f*d1*d2 << "\n\n";
707  errorFlag = -1000;
708  }
709  art::scalarMultiplyDataField<double>(out_c_f_p_d_d, data_c_1, in_f_p_d_d, true);
710  if (std::abs(rst::vectorNorm(out_c_f_p_d_d, NORM_ONE) - (1.0/data_c_1(0,0))*in_f_p_d_d(0,0,0,0)*c*p*f*d1*d2)/rst::vectorNorm(out_c_f_p_d_d, NORM_ONE) > zero) {
711  *outStream << "\n\nINCORRECT scalarMultiplyDataField (6): check result: "
712  << rst::vectorNorm(out_c_f_p_d_d, NORM_ONE) << " != "
713  << (1.0/data_c_1(0,0))*in_f_p_d_d(0,0,0,0)*c*p*f*d1*d2 << "\n\n";
714  errorFlag = -1000;
715  }
716  } // end scope
717 
718  { // start scope
719  *outStream << "\n************ Checking scalarMultiplyDataData, reciprocal=false (branch 1) ************\n";
720 
721  int c=5, p=9, d1=7, d2=13;
722 
723  Kokkos::View<double**> in_c_p("in_c_p", c, p);
724  Kokkos::View<double***> in_c_p_d("in_c_p_d", c, p, d1);
725  Kokkos::View<double****> in_c_p_d_d("in_c_p_d_d", c, p, d1, d2);
726  Kokkos::View<double**> data_c_p("data_c_p", c, p);
727  Kokkos::View<double**> datainv_c_p("datainv_c_p", c, p);
728  Kokkos::View<double**> data_c_1("data_c_1", c, 1);
729  Kokkos::View<double**> datainv_c_1("datainv_c_1", c, 1);
730  Kokkos::View<double**> out_c_p("out_c_p", c, p);
731  Kokkos::View<double**> outi_c_p("outi_c_p", c, p);
732  Kokkos::View<double***> out_c_p_d("out_c_p_d", c, p, d1);
733  Kokkos::View<double***> outi_c_p_d("outi_c_p_d", c, p, d1);
734  Kokkos::View<double****> out_c_p_d_d("out_c_p_d_d", c, p, d1, d2);
735  Kokkos::View<double****> outi_c_p_d_d("outi_c_p_d_d", c, p, d1, d2);
736  double zero = INTREPID_TOL*10000.0;
737 
738  // fill with random numbers
739  for (unsigned int i=0; i<in_c_p.dimension(0); i++)
740  for (unsigned int j=0; j<in_c_p.dimension(1); j++)
741  in_c_p(i,j) = Teuchos::ScalarTraits<double>::random();
742 
743  for (unsigned int i=0; i<in_c_p_d.dimension(0); i++)
744  for (unsigned int j=0; j<in_c_p_d.dimension(1); j++)
745  for (unsigned int k=0; k<in_c_p_d.dimension(2); k++)
746  in_c_p_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
747 
748  for (unsigned int i=0; i<in_c_p_d_d.dimension(0); i++)
749  for (unsigned int j=0; j<in_c_p_d_d.dimension(1); j++)
750  for (unsigned int k=0; k<in_c_p_d_d.dimension(2); k++)
751  for (unsigned int l=0; l<in_c_p_d_d.dimension(3); l++)
752  in_c_p_d_d(i,j,k,l) = Teuchos::ScalarTraits<double>::random();
753 
754  for (unsigned int i=0; i<data_c_p.dimension(0); i++)
755  for (unsigned int j=0; j<data_c_p.dimension(1); j++){
756  data_c_p(i,j) = Teuchos::ScalarTraits<double>::random();
757  datainv_c_p(i,j) = 1.0 / data_c_p(i,j);
758  }
759 
760  for (unsigned int i=0; i<data_c_1.dimension(0); i++)
761  for (unsigned int j=0; j<data_c_1.dimension(1); j++){
762  data_c_1(i,j) = Teuchos::ScalarTraits<double>::random();
763  datainv_c_1(i,j) = 1.0 / data_c_1(i,j);
764  }
765 
766  art::scalarMultiplyDataData<double>(out_c_p, data_c_p, in_c_p);
767  art::scalarMultiplyDataData<double>(outi_c_p, datainv_c_p, out_c_p);
768  rst::subtract(outi_c_p, in_c_p);
769  if (rst::vectorNorm(outi_c_p, NORM_ONE) > zero) {
770  *outStream << "\n\nINCORRECT scalarMultiplyDataData (1): check scalar inverse property\n\n";
771  errorFlag = -1000;
772  }
773  art::scalarMultiplyDataData<double>(out_c_p_d, data_c_p, in_c_p_d);
774  art::scalarMultiplyDataData<double>(outi_c_p_d, datainv_c_p, out_c_p_d);
775  rst::subtract(outi_c_p_d, in_c_p_d);
776  if (rst::vectorNorm(outi_c_p_d, NORM_ONE) > zero) {
777  *outStream << "\n\nINCORRECT scalarMultiplyDataData (2): check scalar inverse property\n\n";
778  errorFlag = -1000;
779  }
780  art::scalarMultiplyDataData<double>(out_c_p_d_d, data_c_p, in_c_p_d_d);
781  art::scalarMultiplyDataData<double>(outi_c_p_d_d, datainv_c_p, out_c_p_d_d);
782  rst::subtract(outi_c_p_d_d, in_c_p_d_d);
783  if (rst::vectorNorm(outi_c_p_d_d, NORM_ONE) > zero) {
784  *outStream << "\n\nINCORRECT scalarMultiplyDataData (3): check scalar inverse property\n\n";
785  errorFlag = -1000;
786  }
787  art::scalarMultiplyDataData<double>(out_c_p_d_d, data_c_1, in_c_p_d_d);
788  art::scalarMultiplyDataData<double>(outi_c_p_d_d, datainv_c_1, out_c_p_d_d);
789  rst::subtract(outi_c_p_d_d, in_c_p_d_d);
790  if (rst::vectorNorm(outi_c_p_d_d, NORM_ONE) > zero) {
791  *outStream << "\n\nINCORRECT scalarMultiplyDataData (4): check scalar inverse property\n\n";
792  errorFlag = -1000;
793  }
794 
795  // fill with constants
796  for (unsigned int i=0; i<in_c_p_d_d.dimension(0); i++)
797  for (unsigned int j=0; j<in_c_p_d_d.dimension(1); j++)
798  for (unsigned int k=0; k<in_c_p_d_d.dimension(2); k++)
799  for (unsigned int l=0; l<in_c_p_d_d.dimension(3); l++)
800  in_c_p_d_d(i,j,k,l) = 1.0;
801 
802 
803  for (unsigned int i=0; i<data_c_p.dimension(0); i++)
804  for (unsigned int j=0; j<data_c_p.dimension(1); j++){
805  data_c_p(i,j) = 5.0;
806  }
807 
808  for (unsigned int i=0; i<data_c_1.dimension(0); i++)
809  for (unsigned int j=0; j<data_c_1.dimension(1); j++){
810  data_c_1(i,j) = 5.0;
811  }
812 
813  art::scalarMultiplyDataData<double>(out_c_p_d_d, data_c_p, in_c_p_d_d);
814  if (std::abs(rst::vectorNorm(out_c_p_d_d, NORM_ONE) - data_c_p(0,0)*in_c_p_d_d(0,0,0,0)*c*p*d1*d2) > zero) {
815  *outStream << "\n\nINCORRECT scalarMultiplyDataData (5): check result: "
816  << rst::vectorNorm(out_c_p_d_d, NORM_ONE) << " != "
817  << data_c_p(0,0)*in_c_p_d_d(0,0,0,0)*c*p*d1*d2 << "\n\n";
818  errorFlag = -1000;
819  }
820  art::scalarMultiplyDataData<double>(out_c_p_d_d, data_c_1, in_c_p_d_d);
821  if (std::abs(rst::vectorNorm(out_c_p_d_d, NORM_ONE) - data_c_1(0,0)*in_c_p_d_d(0,0,0,0)*c*p*d1*d2) > zero) {
822  *outStream << "\n\nINCORRECT scalarMultiplyDataData (6): check result: "
823  << rst::vectorNorm(out_c_p_d_d, NORM_ONE) << " != "
824  << data_c_p(0,0)*in_c_p_d_d(0,0,0,0)*c*p*d1*d2 << "\n\n";
825  errorFlag = -1000;
826  }
827  } // end scope
828 
829  { // start scope
830 
831  // *outStream << "\n************ Checking scalarMultiplyDataData, reciprocal=false (branch 2) ************\n";
832 
833 
834  int c=5, p=9, d1=7, d2=13;
835 
836  Kokkos::View<double*> in_p("in_p", p);
837  Kokkos::View<double**> in_p_d("in_p_d", p, d1);
838  Kokkos::View<double***> in_p_d_d("in_p_d_d", p, d1, d2);
839  Kokkos::View<double**> in_c_p("in_c_p", c, p);
840  Kokkos::View<double***> in_c_p_d("in_c_p_d", c, p, d1);
841  Kokkos::View<double****> in_c_p_d_d("in_c_p_d_d", c, p, d1, d2);
842  Kokkos::View<double**> data_c_p("data_c_p", c, p);
843  Kokkos::View<double**> datainv_c_p("datainv_c_p", c, p);
844  Kokkos::View<double**> data_c_1("data_c_1", c, 1);
845  Kokkos::View<double**> datainv_c_1("datainv_c_1", c, 1);
846  Kokkos::View<double**> data_c_p_one("data_c_p_one", c, p);
847  Kokkos::View<double**> data_c_1_one("data_c_1_one", c, 1);
848  Kokkos::View<double**> out_c_p("out_c_p", c, p);
849  Kokkos::View<double**> outi_c_p("outi_c_p", c, p);
850  Kokkos::View<double***> out_c_p_d("out_c_p_d", c, p, d1);
851  Kokkos::View<double***> outi_c_p_d("outi_c_p_d", c, p, d1);
852  Kokkos::View<double****> out_c_p_d_d("out_c_p_d_d", c, p, d1, d2);
853  Kokkos::View<double****> outi_c_p_d_d("outi_c_p_d_d", c, p, d1, d2);
854  double zero = INTREPID_TOL*10000.0;
855 
856  // fill with random numbers
857  for (unsigned int i=0; i<in_p.dimension(0); i++) {
858  in_p(i) = Teuchos::ScalarTraits<double>::random();
859  }
860 
861  for (unsigned int i=0; i<in_p_d.dimension(0); i++)
862  for (unsigned int j=0; j<in_p_d.dimension(1); j++){
863  in_p_d(i,j) = Teuchos::ScalarTraits<double>::random();
864  }
865 
866  for (unsigned int i=0; i<in_p_d_d.dimension(0); i++)
867  for (unsigned int j=0; j<in_p_d_d.dimension(1); j++)
868  for (unsigned int k=0; k<in_p_d_d.dimension(2); k++)
869  in_p_d_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
870 
871  for (unsigned int i=0; i<data_c_p.dimension(0); i++)
872  for (unsigned int j=0; j<data_c_p.dimension(1); j++){
873  data_c_p(i,j) = Teuchos::ScalarTraits<double>::random();
874  datainv_c_p(i,j) = 1.0 / data_c_p(i,j);
875  data_c_p_one(i,j) = 1.0;
876  }
877 
878  for (unsigned int i=0; i<data_c_1.dimension(0); i++)
879  for (unsigned int j=0; j<data_c_1.dimension(1); j++){
880  data_c_1(i,j) = Teuchos::ScalarTraits<double>::random();
881  datainv_c_1(i,j) = 1.0 / data_c_1(i,j);
882  data_c_1_one(i,j) = 1.0;
883  }
884 
885 
886 
887  art::scalarMultiplyDataData<double>(out_c_p, data_c_p, in_p);
888  art::scalarMultiplyDataData<double>(outi_c_p, datainv_c_p, out_c_p);
889  art::scalarMultiplyDataData<double>(in_c_p, data_c_p_one, in_p);
890  rst::subtract(outi_c_p, in_c_p);
891  if (rst::vectorNorm(outi_c_p, NORM_ONE) > zero) {
892  *outStream << "\n\nINCORRECT scalarMultiplyDataData (1): check scalar inverse property\n\n";
893  errorFlag = -1000;
894  }
895  art::scalarMultiplyDataData<double>(out_c_p_d, data_c_p, in_p_d);
896  art::scalarMultiplyDataData<double>(outi_c_p_d, datainv_c_p, out_c_p_d);
897  art::scalarMultiplyDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
898  rst::subtract(outi_c_p_d, in_c_p_d);
899  if (rst::vectorNorm(outi_c_p_d, NORM_ONE) > zero) {
900  *outStream << "\n\nINCORRECT scalarMultiplyDataData (2): check scalar inverse property\n\n";
901  errorFlag = -1000;
902  }
903  art::scalarMultiplyDataData<double>(out_c_p_d_d, data_c_p, in_p_d_d);
904  art::scalarMultiplyDataData<double>(outi_c_p_d_d, datainv_c_p, out_c_p_d_d);
905  art::scalarMultiplyDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
906  rst::subtract(outi_c_p_d_d, in_c_p_d_d);
907  if (rst::vectorNorm(outi_c_p_d_d, NORM_ONE) > zero) {
908  *outStream << "\n\nINCORRECT scalarMultiplyDataData (3): check scalar inverse property\n\n";
909  errorFlag = -1000;
910  }
911  art::scalarMultiplyDataData<double>(out_c_p_d_d, data_c_1, in_p_d_d);
912  art::scalarMultiplyDataData<double>(outi_c_p_d_d, datainv_c_1, out_c_p_d_d);
913  art::scalarMultiplyDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
914  rst::subtract(outi_c_p_d_d, in_c_p_d_d);
915  if (rst::vectorNorm(outi_c_p_d_d, NORM_ONE) > zero) {
916  *outStream << "\n\nINCORRECT scalarMultiplyDataData (4): check scalar inverse property\n\n";
917  errorFlag = -1000;
918  }
919 
920  // fill with constants
921  for (unsigned int i=0; i<in_p_d_d.dimension(0); i++)
922  for (unsigned int j=0; j<in_p_d_d.dimension(1); j++)
923  for (unsigned int k=0; k<in_p_d_d.dimension(2); k++)
924  in_p_d_d(i,j,k) = 1.0;
925 
926  for (unsigned int i=0; i<data_c_p.dimension(0); i++)
927  for (unsigned int j=0; j<data_c_p.dimension(1); j++){
928  data_c_p(i,j) = 5.0;
929  }
930 
931  for (unsigned int i=0; i<data_c_1.dimension(0); i++)
932  for (unsigned int j=0; j<data_c_1.dimension(1); j++){
933  data_c_1(i,j) = 5.0;
934  }
935 
936  art::scalarMultiplyDataData<double>(out_c_p_d_d, data_c_p, in_p_d_d);
937  if (std::abs(rst::vectorNorm(out_c_p_d_d, NORM_ONE) - data_c_p(0,0)*in_p_d_d(0,0,0)*c*p*d1*d2) > zero) {
938  *outStream << "\n\nINCORRECT scalarMultiplyDataData (5): check result: "
939  << rst::vectorNorm(out_c_p_d_d, NORM_ONE) << " != "
940  << data_c_p(0,0)*in_p_d_d(0,0,0)*c*p*d1*d2 << "\n\n";
941  errorFlag = -1000;
942  }
943  art::scalarMultiplyDataData<double>(out_c_p_d_d, data_c_1, in_p_d_d);
944  if (std::abs(rst::vectorNorm(out_c_p_d_d, NORM_ONE) - data_c_1(0,0)*in_p_d_d(0,0,0)*c*p*d1*d2) > zero) {
945  *outStream << "\n\nINCORRECT scalarMultiplyDataData (6): check result: "
946  << rst::vectorNorm(out_c_p_d_d, NORM_ONE) << " != "
947  << data_c_1(0,0)*in_p_d_d(0,0,0)*c*p*d1*d2 << "\n\n";
948  errorFlag = -1000;
949  }
950  } // end scope
951 
952  { // start scope
953  *outStream << "\n************ Checking scalarMultiplyDataData, reciprocal=true, i.e. division (branch 1) ************\n";
954 
955  int c=5, p=9, d1=7, d2=13;
956 
957  Kokkos::View<double**> in_c_p("in_c_p", c, p);
958  Kokkos::View<double***> in_c_p_d("in_c_p_d", c, p, d1);
959  Kokkos::View<double****> in_c_p_d_d("in_c_p_d_d", c, p, d1, d2);
960  Kokkos::View<double**> data_c_p("data_c_p", c, p);
961  Kokkos::View<double**> datainv_c_p("datainv_c_p", c, p);
962  Kokkos::View<double**> data_c_1("data_c_1", c, 1);
963  Kokkos::View<double**> datainv_c_1("datainv_c_1", c, 1);
964  Kokkos::View<double**> out_c_p("out_c_p", c, p);
965  Kokkos::View<double**> outi_c_p("outi_c_p", c, p);
966  Kokkos::View<double***> out_c_p_d("out_c_p_d", c, p, d1);
967  Kokkos::View<double***> outi_c_p_d("outi_c_p_d", c, p, d1);
968  Kokkos::View<double****> out_c_p_d_d("out_c_p_d_d", c, p, d1, d2);
969  Kokkos::View<double****> outi_c_p_d_d("outi_c_p_d_d", c, p, d1, d2);
970  double zero = INTREPID_TOL*10000.0;
971 
972  // fill with random numbers
973 
974  for (unsigned int i=0; i<in_c_p.dimension(0); i++)
975  for (unsigned int j=0; j<in_c_p.dimension(1); j++){
976  in_c_p(i,j) = Teuchos::ScalarTraits<double>::random();
977  }
978 
979  for (unsigned int i=0; i<in_c_p_d.dimension(0); i++)
980  for (unsigned int j=0; j<in_c_p_d.dimension(1); j++)
981  for (unsigned int k=0; k<in_c_p_d.dimension(2); k++)
982  in_c_p_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
983 
984  for (unsigned int i=0; i<in_c_p_d_d.dimension(0); i++)
985  for (unsigned int j=0; j<in_c_p_d_d.dimension(1); j++)
986  for (unsigned int k=0; k<in_c_p_d_d.dimension(2); k++)
987  for (unsigned int l=0; l<in_c_p_d_d.dimension(3); l++)
988  in_c_p_d_d(i,j,k,l) = Teuchos::ScalarTraits<double>::random();
989 
990  for (unsigned int i=0; i<data_c_p.dimension(0); i++)
991  for (unsigned int j=0; j<data_c_p.dimension(1); j++){
992  data_c_p(i,j) = Teuchos::ScalarTraits<double>::random();
993  datainv_c_p(i,j) = 1.0 / data_c_p(i,j);
994  }
995 
996  for (unsigned int i=0; i<data_c_1.dimension(0); i++)
997  for (unsigned int j=0; j<data_c_1.dimension(1); j++){
998  data_c_1(i,j) = Teuchos::ScalarTraits<double>::random();
999  datainv_c_1(i,j) = 1.0 / data_c_1(i,j);
1000  }
1001 
1002 
1003  art::scalarMultiplyDataData<double>(out_c_p, data_c_p, in_c_p, true);
1004  art::scalarMultiplyDataData<double>(outi_c_p, datainv_c_p, out_c_p, true);
1005  rst::subtract(outi_c_p, in_c_p);
1006  if (rst::vectorNorm(outi_c_p, NORM_ONE) > zero) {
1007  *outStream << "\n\nINCORRECT scalarMultiplyDataData (1): check scalar inverse property\n\n";
1008  errorFlag = -1000;
1009  }
1010  art::scalarMultiplyDataData<double>(out_c_p_d, data_c_p, in_c_p_d, true);
1011  art::scalarMultiplyDataData<double>(outi_c_p_d, datainv_c_p, out_c_p_d, true);
1012  rst::subtract(outi_c_p_d, in_c_p_d);
1013  if (rst::vectorNorm(outi_c_p_d, NORM_ONE) > zero) {
1014  *outStream << "\n\nINCORRECT scalarMultiplyDataData (2): check scalar inverse property\n\n";
1015  errorFlag = -1000;
1016  }
1017  art::scalarMultiplyDataData<double>(out_c_p_d_d, data_c_p, in_c_p_d_d, true);
1018  art::scalarMultiplyDataData<double>(outi_c_p_d_d, datainv_c_p, out_c_p_d_d, true);
1019  rst::subtract(outi_c_p_d_d, in_c_p_d_d);
1020  if (rst::vectorNorm(outi_c_p_d_d, NORM_ONE) > zero) {
1021  *outStream << "\n\nINCORRECT scalarMultiplyDataData (3): check scalar inverse property\n\n";
1022  errorFlag = -1000;
1023  }
1024  art::scalarMultiplyDataData<double>(out_c_p_d_d, data_c_1, in_c_p_d_d, true);
1025  art::scalarMultiplyDataData<double>(outi_c_p_d_d, datainv_c_1, out_c_p_d_d, true);
1026  rst::subtract(outi_c_p_d_d, in_c_p_d_d);
1027  if (rst::vectorNorm(outi_c_p_d_d, NORM_ONE) > zero) {
1028  *outStream << "\n\nINCORRECT scalarMultiplyDataData (4): check scalar inverse property\n\n";
1029  errorFlag = -1000;
1030  }
1031 
1032  // fill with constants
1033  for (unsigned int i=0; i<in_c_p_d_d.dimension(0); i++)
1034  for (unsigned int j=0; j<in_c_p_d_d.dimension(1); j++)
1035  for (unsigned int k=0; k<in_c_p_d_d.dimension(2); k++)
1036  for (unsigned int l=0; l<in_c_p_d_d.dimension(3); l++)
1037  in_c_p_d_d(i,j,k,l) = 1.0;
1038 
1039  for (unsigned int i=0; i<data_c_p.dimension(0); i++)
1040  for (unsigned int j=0; j<data_c_p.dimension(1); j++){
1041  data_c_p(i,j) = 5.0;
1042  }
1043 
1044  for (unsigned int i=0; i<data_c_1.dimension(0); i++)
1045  for (unsigned int j=0; j<data_c_1.dimension(1); j++){
1046  data_c_1(i,j) = 5.0;
1047  }
1048 
1049  art::scalarMultiplyDataData<double>(out_c_p_d_d, data_c_p, in_c_p_d_d, true);
1050  if (std::abs(rst::vectorNorm(out_c_p_d_d, NORM_ONE) - (1.0/data_c_p(0,0))*in_c_p_d_d(0,0,0,0)*c*p*d1*d2)/rst::vectorNorm(out_c_p_d_d, NORM_ONE) > zero) {
1051  *outStream << "\n\nINCORRECT scalarMultiplyDataData (5): check result: "
1052  << rst::vectorNorm(out_c_p_d_d, NORM_ONE) << " != "
1053  << (1.0/data_c_p(0,0))*in_c_p_d_d(0,0,0,0)*c*p*d1*d2 << "\n\n";
1054  errorFlag = -1000;
1055  }
1056  art::scalarMultiplyDataData<double>(out_c_p_d_d, data_c_1, in_c_p_d_d, true);
1057  if (std::abs(rst::vectorNorm(out_c_p_d_d, NORM_ONE) - (1.0/data_c_1(0,0))*in_c_p_d_d(0,0,0,0)*c*p*d1*d2)/rst::vectorNorm(out_c_p_d_d, NORM_ONE) > zero) {
1058  *outStream << "\n\nINCORRECT scalarMultiplyDataData (6): check result: "
1059  << rst::vectorNorm(out_c_p_d_d, NORM_ONE) << " != "
1060  << (1.0/data_c_p(0,0))*in_c_p_d_d(0,0,0,0)*c*p*d1*d2 << "\n\n";
1061  errorFlag = -1000;
1062  }
1063  } // end scope
1064 
1065  { // start scope
1066  *outStream << "\n************ Checking scalarMultiplyDataData, reciprocal=true, i.e. division (branch 2) ************\n";
1067 
1068  int c=5, p=9, d1=7, d2=13;
1069 
1070  Kokkos::View<double*> in_p("in_p", p);
1071  Kokkos::View<double**> in_p_d("in_p_d", p, d1);
1072  Kokkos::View<double***> in_p_d_d("in_p_d_d", p, d1, d2);
1073  Kokkos::View<double**> in_c_p("in_c_p", c, p);
1074  Kokkos::View<double***> in_c_p_d("in_c_p_d", c, p, d1);
1075  Kokkos::View<double****> in_c_p_d_d("in_c_p_d_d", c, p, d1, d2);
1076  Kokkos::View<double**> data_c_p("data_c_p", c, p);
1077  Kokkos::View<double**> datainv_c_p("datainv_c_p", c, p);
1078  Kokkos::View<double**> data_c_1("data_c_1", c, 1);
1079  Kokkos::View<double**> datainv_c_1("datainv_c_1", c, 1);
1080  Kokkos::View<double**> data_c_p_one("data_c_p_one", c, p);
1081  Kokkos::View<double**> data_c_1_one("data_c_1_one", c, 1);
1082  Kokkos::View<double**> out_c_p("out_c_p", c, p);
1083  Kokkos::View<double**> outi_c_p("outi_c_p", c, p);
1084  Kokkos::View<double***> out_c_p_d("out_c_p_d", c, p, d1);
1085  Kokkos::View<double***> outi_c_p_d("outi_c_p_d", c, p, d1);
1086  Kokkos::View<double****> out_c_p_d_d("out_c_p_d_d", c, p, d1, d2);
1087  Kokkos::View<double****> outi_c_p_d_d("outi_c_p_d_d", c, p, d1, d2);
1088  double zero = INTREPID_TOL*10000.0;
1089 
1090  // fill with random numbers
1091  for (unsigned int i=0; i<in_p.dimension(0); i++)
1092  in_p(i) = Teuchos::ScalarTraits<double>::random();
1093 
1094  for (unsigned int i=0; i<in_p_d.dimension(0); i++)
1095  for (unsigned int j=0; j<in_p_d.dimension(1); j++){
1096  in_p_d(i,j) = Teuchos::ScalarTraits<double>::random();
1097  }
1098 
1099  for (unsigned int i=0; i<in_p_d_d.dimension(0); i++)
1100  for (unsigned int j=0; j<in_p_d_d.dimension(1); j++)
1101  for (unsigned int k=0; k<in_p_d_d.dimension(2); k++)
1102  in_p_d_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
1103 
1104  for (unsigned int i=0; i<data_c_p.dimension(0); i++)
1105  for (unsigned int j=0; j<data_c_p.dimension(1); j++){
1106  data_c_p(i,j) = Teuchos::ScalarTraits<double>::random();
1107  datainv_c_p(i,j) = 1.0 / data_c_p(i,j);
1108  data_c_p_one(i,j) = 1.0;
1109  }
1110 
1111  for (unsigned int i=0; i<data_c_1.dimension(0); i++)
1112  for (unsigned int j=0; j<data_c_1.dimension(1); j++){
1113  data_c_1(i,j) = Teuchos::ScalarTraits<double>::random();
1114  datainv_c_1(i,j) = 1.0 / data_c_1(i,j);
1115  data_c_1_one(i,j) = 1.0;
1116  }
1117 
1118  art::scalarMultiplyDataData<double>(out_c_p, data_c_p, in_p, true);
1119  art::scalarMultiplyDataData<double>(outi_c_p, datainv_c_p, out_c_p, true);
1120  art::scalarMultiplyDataData<double>(in_c_p, data_c_p_one, in_p);
1121  rst::subtract(outi_c_p, in_c_p);
1122  if (rst::vectorNorm(outi_c_p, NORM_ONE) > zero) {
1123  *outStream << "\n\nINCORRECT scalarMultiplyDataData (1): check scalar inverse property\n\n";
1124  errorFlag = -1000;
1125  }
1126  art::scalarMultiplyDataData<double>(out_c_p_d, data_c_p, in_p_d, true);
1127  art::scalarMultiplyDataData<double>(outi_c_p_d, datainv_c_p, out_c_p_d, true);
1128  art::scalarMultiplyDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
1129  rst::subtract(outi_c_p_d, in_c_p_d);
1130  if (rst::vectorNorm(outi_c_p_d, NORM_ONE) > zero) {
1131  *outStream << "\n\nINCORRECT scalarMultiplyDataData (2): check scalar inverse property\n\n";
1132  errorFlag = -1000;
1133  }
1134  art::scalarMultiplyDataData<double>(out_c_p_d_d, data_c_p, in_p_d_d, true);
1135  art::scalarMultiplyDataData<double>(outi_c_p_d_d, datainv_c_p, out_c_p_d_d, true);
1136  art::scalarMultiplyDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
1137  rst::subtract(outi_c_p_d_d, in_c_p_d_d);
1138  if (rst::vectorNorm(outi_c_p_d_d, NORM_ONE) > zero) {
1139  *outStream << "\n\nINCORRECT scalarMultiplyDataData (3): check scalar inverse property\n\n";
1140  errorFlag = -1000;
1141  }
1142  art::scalarMultiplyDataData<double>(out_c_p_d_d, data_c_1, in_p_d_d, true);
1143  art::scalarMultiplyDataData<double>(outi_c_p_d_d, datainv_c_1, out_c_p_d_d, true);
1144  art::scalarMultiplyDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
1145  rst::subtract(outi_c_p_d_d, in_c_p_d_d);
1146  if (rst::vectorNorm(outi_c_p_d_d, NORM_ONE) > zero) {
1147  *outStream << "\n\nINCORRECT scalarMultiplyDataData (4): check scalar inverse property\n\n";
1148  errorFlag = -1000;
1149  }
1150 
1151  // fill with constants
1152  for (unsigned int i=0; i<in_p_d_d.dimension(0); i++)
1153  for (unsigned int j=0; j<in_p_d_d.dimension(1); j++)
1154  for (unsigned int k=0; k<in_p_d_d.dimension(2); k++){
1155  in_p_d_d(i,j,k) = 1.0;
1156  }
1157 
1158  for (unsigned int i=0; i<data_c_p.dimension(0); i++)
1159  for (unsigned int j=0; j<data_c_p.dimension(1); j++){
1160  data_c_p(i,j) = 5.0;
1161  }
1162 
1163  for (unsigned int i=0; i<data_c_1.dimension(0); i++)
1164  for (unsigned int j=0; j<data_c_1.dimension(1); j++){
1165  data_c_1(i,j) = 5.0;
1166  }
1167 
1168  art::scalarMultiplyDataData<double>(out_c_p_d_d, data_c_p, in_p_d_d, true);
1169  if (std::abs(rst::vectorNorm(out_c_p_d_d, NORM_ONE) - (1.0/data_c_p(0,0))*in_p_d_d(0,0,0)*c*p*d1*d2)/rst::vectorNorm(out_c_p_d_d, NORM_ONE) > zero) {
1170  *outStream << "\n\nINCORRECT scalarMultiplyDataData (5): check result: "
1171  << rst::vectorNorm(out_c_p_d_d, NORM_ONE) << " != "
1172  << (1.0/data_c_p(0,0))*in_p_d_d(0,0,0)*c*p*d1*d2 << "\n\n";
1173  errorFlag = -1000;
1174  }
1175  art::scalarMultiplyDataData<double>(out_c_p_d_d, data_c_1, in_p_d_d, true);
1176  if (std::abs(rst::vectorNorm(out_c_p_d_d, NORM_ONE) - (1.0/data_c_1(0,0))*in_p_d_d(0,0,0)*c*p*d1*d2)/rst::vectorNorm(out_c_p_d_d, NORM_ONE) > zero) {
1177  *outStream << "\n\nINCORRECT scalarMultiplyDataData (6): check result: "
1178  << rst::vectorNorm(out_c_p_d_d, NORM_ONE) << " != "
1179  << (1.0/data_c_1(0,0))*in_p_d_d(0,0,0)*c*p*d1*d2 << "\n\n";
1180  errorFlag = -1000;
1181  }
1182  } // end scope
1183  /******************************************/
1184  *outStream << "\n";
1185  }
1186  catch (std::logic_error err) {
1187  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
1188  *outStream << err.what() << '\n';
1189  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
1190  errorFlag = -1000;
1191  };
1192 
1193 
1194  if (errorFlag != 0)
1195  std::cout << "End Result: TEST FAILED\n";
1196  else
1197  std::cout << "End Result: TEST PASSED\n";
1198 
1199  // reset format state of std::cout
1200  std::cout.copyfmt(oldFormatState);
1201  Kokkos::finalize();
1202  return errorFlag;
1203 }
Implementation of basic linear algebra functionality in Euclidean space.
Header file for utility class to provide multidimensional containers.
Header file for utility class to provide array tools, such as tensor contractions, etc.
static void scalarMultiplyDataData(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight, const bool reciprocal=false)
There are two use cases: (1) multiplies a rank-2, 3, or 4 container inputDataRight with dimensions (C...
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D...
Utility class that provides methods for higher-order algebraic manipulation of user-defined arrays...
static void scalarMultiplyDataField(ArrayOutFields &outputFields, const ArrayInData &inputData, const ArrayInFields &inputFields, const bool reciprocal=false)
There are two use cases: (1) multiplies a rank-3, 4, or 5 container inputFields with dimensions (C...