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"
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 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 }
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_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::ArrayTools::scalarMultiplyDataData
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...
Definition: Intrepid_ArrayToolsDefScalar.hpp:501
Intrepid::FieldContainer< double >
Intrepid::ArrayTools::scalarMultiplyDataField
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,...
Definition: Intrepid_ArrayToolsDefScalar.hpp:77
Intrepid_ArrayTools.hpp
Header file for utility class to provide array tools, such as tensor contractions,...
Generated by   doxygen 1.8.16