Intrepid
test_02.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 
57 using namespace std;
58 using namespace Intrepid;
59 
60 #define INTREPID_TEST_COMMAND( S ) \
61 { \
62  try { \
63  S ; \
64  } \
65  catch (const 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 
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 
113  typedef ArrayTools art;
114  typedef RealSpaceTools<double> rst;
115 #ifdef HAVE_INTREPID_DEBUG
116  ArrayTools atools;
117 #endif
118 
119  *outStream \
120  << "\n"
121  << "===============================================================================\n"\
122  << "| TEST 1: exceptions |\n"\
123  << "===============================================================================\n";
124 
125  try{
126 
127 #ifdef HAVE_INTREPID_DEBUG
128  FieldContainer<double> a_2_2(2, 2);
129  FieldContainer<double> a_10_2(10, 2);
130  FieldContainer<double> a_10_3(10, 3);
131  FieldContainer<double> a_10_2_2(10, 2, 2);
132  FieldContainer<double> a_10_2_3(10, 2, 3);
133  FieldContainer<double> a_10_3_2(10, 3, 2);
134  FieldContainer<double> a_9_2_2(9, 2, 2);
135 
136  FieldContainer<double> a_10_2_2_2(10, 2, 2, 2);
137  FieldContainer<double> a_9_2_2_2(9, 2, 2, 2);
138  FieldContainer<double> a_10_3_2_2(10, 3, 2, 2);
139  FieldContainer<double> a_10_2_3_2(10, 2, 3, 2);
140  FieldContainer<double> a_10_2_2_3(10, 2, 2, 3);
141 
142  FieldContainer<double> a_10_2_2_2_2(10, 2, 2, 2, 2);
143  FieldContainer<double> a_9_2_2_2_2(9, 2, 2, 2, 2);
144  FieldContainer<double> a_10_3_2_2_2(10, 3, 2, 2, 2);
145  FieldContainer<double> a_10_2_3_2_2(10, 2, 3, 2, 2);
146  FieldContainer<double> a_10_2_2_3_2(10, 2, 2, 3, 2);
147  FieldContainer<double> a_10_2_2_2_3(10, 2, 2, 2, 3);
148 
149  FieldContainer<double> a_9_2(9, 2);
150  FieldContainer<double> a_10_1(10, 1);
151 
152  FieldContainer<double> a_10_1_2(10, 1, 2);
153  FieldContainer<double> a_10_1_3(10, 1, 3);
154 
155  FieldContainer<double> a_10_1_2_2(10, 1, 2, 2);
156 
157  FieldContainer<double> a_2_3_2_2(2, 3, 2, 2);
158  FieldContainer<double> a_2_2_2_2(2, 2, 2, 2);
159  FieldContainer<double> a_2_10(2, 10);
160  FieldContainer<double> a_2(2);
161 
162  *outStream << "-> scalarMultiplyDataField:\n";
163  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_2_2, a_10_2_2, a_2_2) );
164  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_2_2, a_2_2, a_2_2) );
165  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_2_2, a_2_2, a_10_2_2) );
166  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_2_2, a_2_2, a_10_2_2) );
167  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_2_2, a_10_3, a_10_2_2) );
168  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_9_2_2_2_2, a_10_2, a_10_2_2_2_2) );
169  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_3_2_2_2, a_10_2, a_10_2_2_2_2) );
170  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_2_3_2_2, a_10_2, a_10_2_2_2_2) );
171  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_2_2_3_2, a_10_2, a_10_2_2_2_2) );
172  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_2_2_2_3, a_10_2, a_10_2_2_2_2) );
173  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_2_2_2_2, a_10_2, a_10_2_2_2_2) );
174  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_2_2_2_2, a_10_1, a_10_2_2_2_2) );
175  //
176  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_2_2, a_10_2_2, a_2) );
177  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_2_2_2, a_2_2, a_2) );
178  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_2_2_2, a_2_2, a_10_2) );
179  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_2_2, a_10_2, a_2_10) );
180  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_2_2_2_2, a_9_2, a_2_2_2_2) );
181  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_3_2_2_2, a_10_2, a_2_2_2_2) );
182  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_2_3_2_2, a_10_2, a_2_2_2_2) );
183  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_2_2_3_2, a_10_2, a_2_2_2_2) );
184  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_2_2_2_3, a_10_2, a_2_2_2_2) );
185  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_2_2_2_2, a_10_2, a_2_2_2_2) );
186  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataField<double>(a_10_2_2_2_2, a_10_1, a_2_2_2_2) );
187 
188 
189  FieldContainer<double> a_2_2_2(2, 2, 2);
190 
191  *outStream << "-> scalarMultiplyDataData:\n";
192  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_2_2, a_10_2_2, a_2_2) );
193  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_2, a_2_2, a_2) );
194  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_2_2, a_2_2, a_10_2_2) );
195  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_10_2_2, a_2_2, a_10_2_2) );
196  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_10_2_2, a_10_3, a_10_2_2) );
197  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_9_2_2_2, a_10_2, a_10_2_2_2) );
198  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_10_3_2_2, a_10_2, a_10_2_2_2) );
199  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_10_2_3_2, a_10_2, a_10_2_2_2) );
200  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_10_2_2_3, a_10_2, a_10_2_2_2) );
201  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_10_2_2_2, a_10_2, a_10_2_2_2) );
202  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_10_2_2_2, a_10_1, a_10_2_2_2) );
203  //
204  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_2_2, a_10_2_2, a_2) );
205  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_10_2_2_2_2, a_2_2, a_10_2_2_2) );
206  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_10_2_2_2, a_2_2, a_10_2) );
207  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_10_2_2, a_2_2, a_10_2) );
208  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_10_2_2_2, a_9_2, a_2_2_2) );
209  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_10_3_2_2, a_10_2, a_2_2_2) );
210  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_10_2_3_2, a_10_2, a_2_2_2) );
211  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_10_2_2_3, a_10_2, a_2_2_2) );
212  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_10_2_2_2, a_10_2, a_2_2_2) );
213  INTREPID_TEST_COMMAND( atools.scalarMultiplyDataData<double>(a_10_2_2_2, a_10_1, a_2_2_2) );
214 #endif
215 
216  }
217  catch (const std::logic_error & err) {
218  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
219  *outStream << err.what() << '\n';
220  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
221  errorFlag = -1000;
222  };
223 
224 #ifdef HAVE_INTREPID_DEBUG
225  if (Teuchos::TestForException_getThrowNumber() != endThrowNumber)
226  errorFlag++;
227 #endif
228 
229  *outStream \
230  << "\n"
231  << "===============================================================================\n"\
232  << "| TEST 2: correctness of math operations |\n"\
233  << "===============================================================================\n";
234 
235  outStream->precision(20);
236 
237  try {
238  { // start scope
239  *outStream << "\n************ Checking scalarMultiplyDataField, reciprocal=false (branch 1) ************\n";
240 
241  int c=5, p=9, f=7, d1=7, d2=13;
242 
243  FieldContainer<double> in_c_f_p(c, f, p);
244  FieldContainer<double> in_c_f_p_d(c, f, p, d1);
245  FieldContainer<double> in_c_f_p_d_d(c, f, p, d1, d2);
246  FieldContainer<double> data_c_p(c, p);
247  FieldContainer<double> datainv_c_p(c, p);
248  FieldContainer<double> data_c_1(c, 1);
249  FieldContainer<double> datainv_c_1(c, 1);
250  FieldContainer<double> out_c_f_p(c, f, p);
251  FieldContainer<double> outi_c_f_p(c, f, p);
252  FieldContainer<double> out_c_f_p_d(c, f, p, d1);
253  FieldContainer<double> outi_c_f_p_d(c, f, p, d1);
254  FieldContainer<double> out_c_f_p_d_d(c, f, p, d1, d2);
255  FieldContainer<double> outi_c_f_p_d_d(c, f, p, d1, d2);
256  double zero = INTREPID_TOL*10000.0;
257 
258  // fill with random numbers
259  for (int i=0; i<in_c_f_p.size(); i++) {
260  in_c_f_p[i] = Teuchos::ScalarTraits<double>::random();
261  }
262  for (int i=0; i<in_c_f_p_d.size(); i++) {
263  in_c_f_p_d[i] = Teuchos::ScalarTraits<double>::random();
264  }
265  for (int i=0; i<in_c_f_p_d_d.size(); i++) {
266  in_c_f_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
267  }
268  for (int i=0; i<data_c_p.size(); i++) {
269  data_c_p[i] = Teuchos::ScalarTraits<double>::random();
270  datainv_c_p[i] = 1.0 / data_c_p[i];
271  }
272  for (int i=0; i<data_c_1.size(); i++) {
273  data_c_1[i] = Teuchos::ScalarTraits<double>::random();
274  datainv_c_1[i] = 1.0 / data_c_1[i];
275  }
276 
277  art::scalarMultiplyDataField<double>(out_c_f_p, data_c_p, in_c_f_p);
278  art::scalarMultiplyDataField<double>(outi_c_f_p, datainv_c_p, out_c_f_p);
279  rst::subtract(&outi_c_f_p[0], &in_c_f_p[0], outi_c_f_p.size());
280  if (rst::vectorNorm(&outi_c_f_p[0], outi_c_f_p.size(), NORM_ONE) > zero) {
281  *outStream << "\n\nINCORRECT scalarMultiplyDataField (1): check scalar inverse property\n\n";
282  errorFlag = -1000;
283  }
284  art::scalarMultiplyDataField<double>(out_c_f_p_d, data_c_p, in_c_f_p_d);
285  art::scalarMultiplyDataField<double>(outi_c_f_p_d, datainv_c_p, out_c_f_p_d);
286  rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
287  if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
288  *outStream << "\n\nINCORRECT scalarMultiplyDataField (2): check scalar inverse property\n\n";
289  errorFlag = -1000;
290  }
291  art::scalarMultiplyDataField<double>(out_c_f_p_d_d, data_c_p, in_c_f_p_d_d);
292  art::scalarMultiplyDataField<double>(outi_c_f_p_d_d, datainv_c_p, out_c_f_p_d_d);
293  rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
294  if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
295  *outStream << "\n\nINCORRECT scalarMultiplyDataField (3): check scalar inverse property\n\n";
296  errorFlag = -1000;
297  }
298  art::scalarMultiplyDataField<double>(out_c_f_p_d_d, data_c_1, in_c_f_p_d_d);
299  art::scalarMultiplyDataField<double>(outi_c_f_p_d_d, datainv_c_1, out_c_f_p_d_d);
300  rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
301  if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
302  *outStream << "\n\nINCORRECT scalarMultiplyDataField (4): check scalar inverse property\n\n";
303  errorFlag = -1000;
304  }
305 
306  // fill with constants
307  for (int i=0; i<in_c_f_p_d_d.size(); i++) {
308  in_c_f_p_d_d[i] = 1.0;
309  }
310  for (int i=0; i<data_c_p.size(); i++) {
311  data_c_p[i] = 5.0;
312  }
313  for (int i=0; i<data_c_1.size(); i++) {
314  data_c_1[i] = 5.0;
315  }
316  art::scalarMultiplyDataField<double>(out_c_f_p_d_d, data_c_p, in_c_f_p_d_d);
317  if (std::abs(rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) - data_c_p[0]*in_c_f_p_d_d[0]*c*p*f*d1*d2) > zero) {
318  *outStream << "\n\nINCORRECT scalarMultiplyDataField (5): check result: "
319  << rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) << " != "
320  << data_c_p[0]*in_c_f_p_d_d[0]*c*f*p*d1*d2 << "\n\n";
321  errorFlag = -1000;
322  }
323  art::scalarMultiplyDataField<double>(out_c_f_p_d_d, data_c_1, in_c_f_p_d_d);
324  if (std::abs(rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) - data_c_1[0]*in_c_f_p_d_d[0]*c*p*f*d1*d2) > zero) {
325  *outStream << "\n\nINCORRECT scalarMultiplyDataField (6): check result: "
326  << rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) << " != "
327  << data_c_p[0]*in_c_f_p_d_d[0]*c*f*p*d1*d2 << "\n\n";
328  errorFlag = -1000;
329  }
330  } // end scope
331 
332  { // start scope
333  *outStream << "\n************ Checking scalarMultiplyDataField, reciprocal=false (branch 2) ************\n";
334 
335  int c=5, p=9, f=7, d1=7, d2=13;
336 
337  FieldContainer<double> in_f_p(f, p);
338  FieldContainer<double> in_f_p_d(f, p, d1);
339  FieldContainer<double> in_f_p_d_d(f, p, d1, d2);
340  FieldContainer<double> in_c_f_p(c, f, p);
341  FieldContainer<double> in_c_f_p_d(c, f, p, d1);
342  FieldContainer<double> in_c_f_p_d_d(c, f, p, d1, d2);
343  FieldContainer<double> data_c_p(c, p);
344  FieldContainer<double> datainv_c_p(c, p);
345  FieldContainer<double> data_c_1(c, 1);
346  FieldContainer<double> datainv_c_1(c, 1);
347  FieldContainer<double> data_c_p_one(c, p);
348  FieldContainer<double> data_c_1_one(c, 1);
349  FieldContainer<double> out_c_f_p(c, f, p);
350  FieldContainer<double> outi_c_f_p(c, f, p);
351  FieldContainer<double> out_c_f_p_d(c, f, p, d1);
352  FieldContainer<double> outi_c_f_p_d(c, f, p, d1);
353  FieldContainer<double> out_c_f_p_d_d(c, f, p, d1, d2);
354  FieldContainer<double> outi_c_f_p_d_d(c, f, p, d1, d2);
355  double zero = INTREPID_TOL*10000.0;
356 
357  // fill with random numbers
358  for (int i=0; i<in_f_p.size(); i++) {
359  in_f_p[i] = Teuchos::ScalarTraits<double>::random();
360  }
361  for (int i=0; i<in_f_p_d.size(); i++) {
362  in_f_p_d[i] = Teuchos::ScalarTraits<double>::random();
363  }
364  for (int i=0; i<in_f_p_d_d.size(); i++) {
365  in_f_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
366  }
367  for (int i=0; i<data_c_p.size(); i++) {
368  data_c_p[i] = Teuchos::ScalarTraits<double>::random();
369  datainv_c_p[i] = 1.0 / data_c_p[i];
370  data_c_p_one[i] = 1.0;
371  }
372  for (int i=0; i<data_c_1.size(); i++) {
373  data_c_1[i] = Teuchos::ScalarTraits<double>::random();
374  datainv_c_1[i] = 1.0 / data_c_1[i];
375  data_c_1_one[i] = 1.0;
376  }
377 
378  art::scalarMultiplyDataField<double>(out_c_f_p, data_c_p, in_f_p);
379  art::scalarMultiplyDataField<double>(outi_c_f_p, datainv_c_p, out_c_f_p);
380  art::scalarMultiplyDataField<double>(in_c_f_p, data_c_p_one, in_f_p);
381  rst::subtract(&outi_c_f_p[0], &in_c_f_p[0], outi_c_f_p.size());
382  if (rst::vectorNorm(&outi_c_f_p[0], outi_c_f_p.size(), NORM_ONE) > zero) {
383  *outStream << "\n\nINCORRECT scalarMultiplyDataField (1): check scalar inverse property\n\n";
384  errorFlag = -1000;
385  }
386  art::scalarMultiplyDataField<double>(out_c_f_p_d, data_c_p, in_f_p_d);
387  art::scalarMultiplyDataField<double>(outi_c_f_p_d, datainv_c_p, out_c_f_p_d);
388  art::scalarMultiplyDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
389  rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
390  if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
391  *outStream << "\n\nINCORRECT scalarMultiplyDataField (2): check scalar inverse property\n\n";
392  errorFlag = -1000;
393  }
394  art::scalarMultiplyDataField<double>(out_c_f_p_d_d, data_c_p, in_f_p_d_d);
395  art::scalarMultiplyDataField<double>(outi_c_f_p_d_d, datainv_c_p, out_c_f_p_d_d);
396  art::scalarMultiplyDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
397  rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
398  if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
399  *outStream << "\n\nINCORRECT scalarMultiplyDataField (3): check scalar inverse property\n\n";
400  errorFlag = -1000;
401  }
402  art::scalarMultiplyDataField<double>(out_c_f_p_d_d, data_c_1, in_f_p_d_d);
403  art::scalarMultiplyDataField<double>(outi_c_f_p_d_d, datainv_c_1, out_c_f_p_d_d);
404  art::scalarMultiplyDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
405  rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
406  if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
407  *outStream << "\n\nINCORRECT scalarMultiplyDataField (4): check scalar inverse property\n\n";
408  errorFlag = -1000;
409  }
410 
411  // fill with constants
412  for (int i=0; i<in_f_p_d_d.size(); i++) {
413  in_f_p_d_d[i] = 1.0;
414  }
415  for (int i=0; i<data_c_p.size(); i++) {
416  data_c_p[i] = 5.0;
417  }
418  for (int i=0; i<data_c_1.size(); i++) {
419  data_c_1[i] = 5.0;
420  }
421  art::scalarMultiplyDataField<double>(out_c_f_p_d_d, data_c_p, in_f_p_d_d);
422  if (std::abs(rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) - data_c_p[0]*in_f_p_d_d[0]*c*p*f*d1*d2) > zero) {
423  *outStream << "\n\nINCORRECT scalarMultiplyDataField (5): check result: "
424  << rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) << " != "
425  << data_c_p[0]*in_f_p_d_d[0]*c*f*p*d1*d2 << "\n\n";
426  errorFlag = -1000;
427  }
428  art::scalarMultiplyDataField<double>(out_c_f_p_d_d, data_c_1, in_f_p_d_d);
429  if (std::abs(rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) - data_c_1[0]*in_f_p_d_d[0]*c*p*f*d1*d2) > zero) {
430  *outStream << "\n\nINCORRECT scalarMultiplyDataField (6): check result: "
431  << rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) << " != "
432  << data_c_1[0]*in_f_p_d_d[0]*c*f*p*d1*d2 << "\n\n";
433  errorFlag = -1000;
434  }
435  } // end scope
436 
437  { // start scope
438  *outStream << "\n************ Checking scalarMultiplyDataField, reciprocal=true, i.e. division (branch 1) ************\n";
439 
440  int c=5, p=9, f=7, d1=7, d2=13;
441 
442  FieldContainer<double> in_c_f_p(c, f, p);
443  FieldContainer<double> in_c_f_p_d(c, f, p, d1);
444  FieldContainer<double> in_c_f_p_d_d(c, f, p, d1, d2);
445  FieldContainer<double> data_c_p(c, p);
446  FieldContainer<double> datainv_c_p(c, p);
447  FieldContainer<double> data_c_1(c, 1);
448  FieldContainer<double> datainv_c_1(c, 1);
449  FieldContainer<double> out_c_f_p(c, f, p);
450  FieldContainer<double> outi_c_f_p(c, f, p);
451  FieldContainer<double> out_c_f_p_d(c, f, p, d1);
452  FieldContainer<double> outi_c_f_p_d(c, f, p, d1);
453  FieldContainer<double> out_c_f_p_d_d(c, f, p, d1, d2);
454  FieldContainer<double> outi_c_f_p_d_d(c, f, p, d1, d2);
455  double zero = INTREPID_TOL*10000.0;
456 
457  // fill with random numbers
458  for (int i=0; i<in_c_f_p.size(); i++) {
459  in_c_f_p[i] = Teuchos::ScalarTraits<double>::random();
460  }
461  for (int i=0; i<in_c_f_p_d.size(); i++) {
462  in_c_f_p_d[i] = Teuchos::ScalarTraits<double>::random();
463  }
464  for (int i=0; i<in_c_f_p_d_d.size(); i++) {
465  in_c_f_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
466  }
467  for (int i=0; i<data_c_p.size(); i++) {
468  data_c_p[i] = Teuchos::ScalarTraits<double>::random();
469  datainv_c_p[i] = 1.0 / data_c_p[i];
470  }
471  for (int i=0; i<data_c_1.size(); i++) {
472  data_c_1[i] = Teuchos::ScalarTraits<double>::random();
473  datainv_c_1[i] = 1.0 / data_c_1[i];
474  }
475 
476  art::scalarMultiplyDataField<double>(out_c_f_p, data_c_p, in_c_f_p, true);
477  art::scalarMultiplyDataField<double>(outi_c_f_p, datainv_c_p, out_c_f_p, true);
478  rst::subtract(&outi_c_f_p[0], &in_c_f_p[0], outi_c_f_p.size());
479  if (rst::vectorNorm(&outi_c_f_p[0], outi_c_f_p.size(), NORM_ONE) > zero) {
480  *outStream << "\n\nINCORRECT scalarMultiplyDataField (1): check scalar inverse property\n\n";
481  errorFlag = -1000;
482  }
483  art::scalarMultiplyDataField<double>(out_c_f_p_d, data_c_p, in_c_f_p_d, true);
484  art::scalarMultiplyDataField<double>(outi_c_f_p_d, datainv_c_p, out_c_f_p_d, true);
485  rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
486  if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
487  *outStream << "\n\nINCORRECT scalarMultiplyDataField (2): check scalar inverse property\n\n";
488  errorFlag = -1000;
489  }
490  art::scalarMultiplyDataField<double>(out_c_f_p_d_d, data_c_p, in_c_f_p_d_d, true);
491  art::scalarMultiplyDataField<double>(outi_c_f_p_d_d, datainv_c_p, out_c_f_p_d_d, true);
492  rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
493  if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
494  *outStream << "\n\nINCORRECT scalarMultiplyDataField (3): check scalar inverse property\n\n";
495  errorFlag = -1000;
496  }
497  art::scalarMultiplyDataField<double>(out_c_f_p_d_d, data_c_1, in_c_f_p_d_d, true);
498  art::scalarMultiplyDataField<double>(outi_c_f_p_d_d, datainv_c_1, out_c_f_p_d_d, true);
499  rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
500  if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
501  *outStream << "\n\nINCORRECT scalarMultiplyDataField (4): check scalar inverse property\n\n";
502  errorFlag = -1000;
503  }
504 
505  // fill with constants
506  for (int i=0; i<in_c_f_p_d_d.size(); i++) {
507  in_c_f_p_d_d[i] = 1.0;
508  }
509  for (int i=0; i<data_c_p.size(); i++) {
510  data_c_p[i] = 5.0;
511  }
512  for (int i=0; i<data_c_1.size(); i++) {
513  data_c_1[i] = 5.0;
514  }
515  art::scalarMultiplyDataField<double>(out_c_f_p_d_d, data_c_p, in_c_f_p_d_d, true);
516  if (std::abs(rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) - (1.0/data_c_p[0])*in_c_f_p_d_d[0]*c*p*f*d1*d2)/rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) > zero) {
517  *outStream << "\n\nINCORRECT scalarMultiplyDataField (5): check result: "
518  << rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) << " != "
519  << (1.0/data_c_p[0])*in_c_f_p_d_d[0]*c*p*f*d1*d2 << "\n\n";
520  errorFlag = -1000;
521  }
522  art::scalarMultiplyDataField<double>(out_c_f_p_d_d, data_c_1, in_c_f_p_d_d, true);
523  if (std::abs(rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) - (1.0/data_c_1[0])*in_c_f_p_d_d[0]*c*p*f*d1*d2)/rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) > zero) {
524  *outStream << "\n\nINCORRECT scalarMultiplyDataField (6): check result: "
525  << rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) << " != "
526  << (1.0/data_c_p[0])*in_c_f_p_d_d[0]*c*p*f*d1*d2 << "\n\n";
527  errorFlag = -1000;
528  }
529  } // end scope
530 
531  { // start scope
532  *outStream << "\n************ Checking scalarMultiplyDataField, reciprocal=true, i.e. division (branch 2) ************\n";
533 
534  int c=5, p=9, f=7, d1=7, d2=13;
535 
536  FieldContainer<double> in_f_p(f, p);
537  FieldContainer<double> in_f_p_d(f, p, d1);
538  FieldContainer<double> in_f_p_d_d(f, p, d1, d2);
539  FieldContainer<double> in_c_f_p(c, f, p);
540  FieldContainer<double> in_c_f_p_d(c, f, p, d1);
541  FieldContainer<double> in_c_f_p_d_d(c, f, p, d1, d2);
542  FieldContainer<double> data_c_p(c, p);
543  FieldContainer<double> datainv_c_p(c, p);
544  FieldContainer<double> data_c_1(c, 1);
545  FieldContainer<double> datainv_c_1(c, 1);
546  FieldContainer<double> data_c_p_one(c, p);
547  FieldContainer<double> data_c_1_one(c, 1);
548  FieldContainer<double> out_c_f_p(c, f, p);
549  FieldContainer<double> outi_c_f_p(c, f, p);
550  FieldContainer<double> out_c_f_p_d(c, f, p, d1);
551  FieldContainer<double> outi_c_f_p_d(c, f, p, d1);
552  FieldContainer<double> out_c_f_p_d_d(c, f, p, d1, d2);
553  FieldContainer<double> outi_c_f_p_d_d(c, f, p, d1, d2);
554  double zero = INTREPID_TOL*10000.0;
555 
556  // fill with random numbers
557  for (int i=0; i<in_f_p.size(); i++) {
558  in_f_p[i] = Teuchos::ScalarTraits<double>::random();
559  }
560  for (int i=0; i<in_f_p_d.size(); i++) {
561  in_f_p_d[i] = Teuchos::ScalarTraits<double>::random();
562  }
563  for (int i=0; i<in_f_p_d_d.size(); i++) {
564  in_f_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
565  }
566  for (int i=0; i<data_c_p.size(); i++) {
567  data_c_p[i] = Teuchos::ScalarTraits<double>::random();
568  datainv_c_p[i] = 1.0 / data_c_p[i];
569  data_c_p_one[i] = 1.0;
570  }
571  for (int i=0; i<data_c_1.size(); i++) {
572  data_c_1[i] = Teuchos::ScalarTraits<double>::random();
573  datainv_c_1[i] = 1.0 / data_c_1[i];
574  data_c_1_one[i] = 1.0;
575  }
576 
577  art::scalarMultiplyDataField<double>(out_c_f_p, data_c_p, in_f_p, true);
578  art::scalarMultiplyDataField<double>(outi_c_f_p, datainv_c_p, out_c_f_p, true);
579  art::scalarMultiplyDataField<double>(in_c_f_p, data_c_p_one, in_f_p);
580  rst::subtract(&outi_c_f_p[0], &in_c_f_p[0], outi_c_f_p.size());
581  if (rst::vectorNorm(&outi_c_f_p[0], outi_c_f_p.size(), NORM_ONE) > zero) {
582  *outStream << "\n\nINCORRECT scalarMultiplyDataField (1): check scalar inverse property\n\n";
583  errorFlag = -1000;
584  }
585  art::scalarMultiplyDataField<double>(out_c_f_p_d, data_c_p, in_f_p_d, true);
586  art::scalarMultiplyDataField<double>(outi_c_f_p_d, datainv_c_p, out_c_f_p_d, true);
587  art::scalarMultiplyDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
588  rst::subtract(&outi_c_f_p_d[0], &in_c_f_p_d[0], outi_c_f_p_d.size());
589  if (rst::vectorNorm(&outi_c_f_p_d[0], outi_c_f_p_d.size(), NORM_ONE) > zero) {
590  *outStream << "\n\nINCORRECT scalarMultiplyDataField (2): check scalar inverse property\n\n";
591  errorFlag = -1000;
592  }
593  art::scalarMultiplyDataField<double>(out_c_f_p_d_d, data_c_p, in_f_p_d_d, true);
594  art::scalarMultiplyDataField<double>(outi_c_f_p_d_d, datainv_c_p, out_c_f_p_d_d, true);
595  art::scalarMultiplyDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
596  rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
597  if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
598  *outStream << "\n\nINCORRECT scalarMultiplyDataField (3): check scalar inverse property\n\n";
599  errorFlag = -1000;
600  }
601  art::scalarMultiplyDataField<double>(out_c_f_p_d_d, data_c_1, in_f_p_d_d, true);
602  art::scalarMultiplyDataField<double>(outi_c_f_p_d_d, datainv_c_1, out_c_f_p_d_d, true);
603  art::scalarMultiplyDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
604  rst::subtract(&outi_c_f_p_d_d[0], &in_c_f_p_d_d[0], outi_c_f_p_d_d.size());
605  if (rst::vectorNorm(&outi_c_f_p_d_d[0], outi_c_f_p_d_d.size(), NORM_ONE) > zero) {
606  *outStream << "\n\nINCORRECT scalarMultiplyDataField (4): check scalar inverse property\n\n";
607  errorFlag = -1000;
608  }
609 
610  // fill with constants
611  for (int i=0; i<in_f_p_d_d.size(); i++) {
612  in_f_p_d_d[i] = 1.0;
613  }
614  for (int i=0; i<data_c_p.size(); i++) {
615  data_c_p[i] = 5.0;
616  }
617  for (int i=0; i<data_c_1.size(); i++) {
618  data_c_1[i] = 5.0;
619  }
620  art::scalarMultiplyDataField<double>(out_c_f_p_d_d, data_c_p, in_f_p_d_d, true);
621  if (std::abs(rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) - (1.0/data_c_p[0])*in_f_p_d_d[0]*c*p*f*d1*d2)/rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) > zero) {
622  *outStream << "\n\nINCORRECT scalarMultiplyDataField (5): check result: "
623  << rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) << " != "
624  << (1.0/data_c_p[0])*in_f_p_d_d[0]*c*p*f*d1*d2 << "\n\n";
625  errorFlag = -1000;
626  }
627  art::scalarMultiplyDataField<double>(out_c_f_p_d_d, data_c_1, in_f_p_d_d, true);
628  if (std::abs(rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) - (1.0/data_c_1[0])*in_f_p_d_d[0]*c*p*f*d1*d2)/rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) > zero) {
629  *outStream << "\n\nINCORRECT scalarMultiplyDataField (6): check result: "
630  << rst::vectorNorm(&out_c_f_p_d_d[0], out_c_f_p_d_d.size(), NORM_ONE) << " != "
631  << (1.0/data_c_1[0])*in_f_p_d_d[0]*c*p*f*d1*d2 << "\n\n";
632  errorFlag = -1000;
633  }
634  } // end scope
635 
636 
637 
638 
639  { // start scope
640  *outStream << "\n************ Checking scalarMultiplyDataData, reciprocal=false (branch 1) ************\n";
641 
642  int c=5, p=9, d1=7, d2=13;
643 
644  FieldContainer<double> in_c_p(c, p);
645  FieldContainer<double> in_c_p_d(c, p, d1);
646  FieldContainer<double> in_c_p_d_d(c, p, d1, d2);
647  FieldContainer<double> data_c_p(c, p);
648  FieldContainer<double> datainv_c_p(c, p);
649  FieldContainer<double> data_c_1(c, 1);
650  FieldContainer<double> datainv_c_1(c, 1);
651  FieldContainer<double> out_c_p(c, p);
652  FieldContainer<double> outi_c_p(c, p);
653  FieldContainer<double> out_c_p_d(c, p, d1);
654  FieldContainer<double> outi_c_p_d(c, p, d1);
655  FieldContainer<double> out_c_p_d_d(c, p, d1, d2);
656  FieldContainer<double> outi_c_p_d_d(c, p, d1, d2);
657  double zero = INTREPID_TOL*10000.0;
658 
659  // fill with random numbers
660  for (int i=0; i<in_c_p.size(); i++) {
661  in_c_p[i] = Teuchos::ScalarTraits<double>::random();
662  }
663  for (int i=0; i<in_c_p_d.size(); i++) {
664  in_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
665  }
666  for (int i=0; i<in_c_p_d_d.size(); i++) {
667  in_c_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
668  }
669  for (int i=0; i<data_c_p.size(); i++) {
670  data_c_p[i] = Teuchos::ScalarTraits<double>::random();
671  datainv_c_p[i] = 1.0 / data_c_p[i];
672  }
673  for (int i=0; i<data_c_1.size(); i++) {
674  data_c_1[i] = Teuchos::ScalarTraits<double>::random();
675  datainv_c_1[i] = 1.0 / data_c_1[i];
676  }
677 
678  art::scalarMultiplyDataData<double>(out_c_p, data_c_p, in_c_p);
679  art::scalarMultiplyDataData<double>(outi_c_p, datainv_c_p, out_c_p);
680  rst::subtract(&outi_c_p[0], &in_c_p[0], outi_c_p.size());
681  if (rst::vectorNorm(&outi_c_p[0], outi_c_p.size(), NORM_ONE) > zero) {
682  *outStream << "\n\nINCORRECT scalarMultiplyDataData (1): check scalar inverse property\n\n";
683  errorFlag = -1000;
684  }
685  art::scalarMultiplyDataData<double>(out_c_p_d, data_c_p, in_c_p_d);
686  art::scalarMultiplyDataData<double>(outi_c_p_d, datainv_c_p, out_c_p_d);
687  rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
688  if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
689  *outStream << "\n\nINCORRECT scalarMultiplyDataData (2): check scalar inverse property\n\n";
690  errorFlag = -1000;
691  }
692  art::scalarMultiplyDataData<double>(out_c_p_d_d, data_c_p, in_c_p_d_d);
693  art::scalarMultiplyDataData<double>(outi_c_p_d_d, datainv_c_p, out_c_p_d_d);
694  rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
695  if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
696  *outStream << "\n\nINCORRECT scalarMultiplyDataData (3): check scalar inverse property\n\n";
697  errorFlag = -1000;
698  }
699  art::scalarMultiplyDataData<double>(out_c_p_d_d, data_c_1, in_c_p_d_d);
700  art::scalarMultiplyDataData<double>(outi_c_p_d_d, datainv_c_1, out_c_p_d_d);
701  rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
702  if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
703  *outStream << "\n\nINCORRECT scalarMultiplyDataData (4): check scalar inverse property\n\n";
704  errorFlag = -1000;
705  }
706 
707  // fill with constants
708  for (int i=0; i<in_c_p_d_d.size(); i++) {
709  in_c_p_d_d[i] = 1.0;
710  }
711  for (int i=0; i<data_c_p.size(); i++) {
712  data_c_p[i] = 5.0;
713  }
714  for (int i=0; i<data_c_1.size(); i++) {
715  data_c_1[i] = 5.0;
716  }
717  art::scalarMultiplyDataData<double>(out_c_p_d_d, data_c_p, in_c_p_d_d);
718  if (std::abs(rst::vectorNorm(&out_c_p_d_d[0], out_c_p_d_d.size(), NORM_ONE) - data_c_p[0]*in_c_p_d_d[0]*c*p*d1*d2) > zero) {
719  *outStream << "\n\nINCORRECT scalarMultiplyDataData (5): check result: "
720  << rst::vectorNorm(&out_c_p_d_d[0], out_c_p_d_d.size(), NORM_ONE) << " != "
721  << data_c_p[0]*in_c_p_d_d[0]*c*p*d1*d2 << "\n\n";
722  errorFlag = -1000;
723  }
724  art::scalarMultiplyDataData<double>(out_c_p_d_d, data_c_1, in_c_p_d_d);
725  if (std::abs(rst::vectorNorm(&out_c_p_d_d[0], out_c_p_d_d.size(), NORM_ONE) - data_c_1[0]*in_c_p_d_d[0]*c*p*d1*d2) > zero) {
726  *outStream << "\n\nINCORRECT scalarMultiplyDataData (6): check result: "
727  << rst::vectorNorm(&out_c_p_d_d[0], out_c_p_d_d.size(), NORM_ONE) << " != "
728  << data_c_p[0]*in_c_p_d_d[0]*c*p*d1*d2 << "\n\n";
729  errorFlag = -1000;
730  }
731  } // end scope
732 
733  { // start scope
734  *outStream << "\n************ Checking scalarMultiplyDataData, reciprocal=false (branch 2) ************\n";
735 
736  int c=5, p=9, d1=7, d2=13;
737 
738  FieldContainer<double> in_p(p);
739  FieldContainer<double> in_p_d(p, d1);
740  FieldContainer<double> in_p_d_d(p, d1, d2);
741  FieldContainer<double> in_c_p(c, p);
742  FieldContainer<double> in_c_p_d(c, p, d1);
743  FieldContainer<double> in_c_p_d_d(c, p, d1, d2);
744  FieldContainer<double> data_c_p(c, p);
745  FieldContainer<double> datainv_c_p(c, p);
746  FieldContainer<double> data_c_1(c, 1);
747  FieldContainer<double> datainv_c_1(c, 1);
748  FieldContainer<double> data_c_p_one(c, p);
749  FieldContainer<double> data_c_1_one(c, 1);
750  FieldContainer<double> out_c_p(c, p);
751  FieldContainer<double> outi_c_p(c, p);
752  FieldContainer<double> out_c_p_d(c, p, d1);
753  FieldContainer<double> outi_c_p_d(c, p, d1);
754  FieldContainer<double> out_c_p_d_d(c, p, d1, d2);
755  FieldContainer<double> outi_c_p_d_d(c, p, d1, d2);
756  double zero = INTREPID_TOL*10000.0;
757 
758  // fill with random numbers
759  for (int i=0; i<in_p.size(); i++) {
760  in_p[i] = Teuchos::ScalarTraits<double>::random();
761  }
762  for (int i=0; i<in_p_d.size(); i++) {
763  in_p_d[i] = Teuchos::ScalarTraits<double>::random();
764  }
765  for (int i=0; i<in_p_d_d.size(); i++) {
766  in_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
767  }
768  for (int i=0; i<data_c_p.size(); i++) {
769  data_c_p[i] = Teuchos::ScalarTraits<double>::random();
770  datainv_c_p[i] = 1.0 / data_c_p[i];
771  data_c_p_one[i] = 1.0;
772  }
773  for (int i=0; i<data_c_1.size(); i++) {
774  data_c_1[i] = Teuchos::ScalarTraits<double>::random();
775  datainv_c_1[i] = 1.0 / data_c_1[i];
776  data_c_1_one[i] = 1.0;
777  }
778 
779  art::scalarMultiplyDataData<double>(out_c_p, data_c_p, in_p);
780  art::scalarMultiplyDataData<double>(outi_c_p, datainv_c_p, out_c_p);
781  art::scalarMultiplyDataData<double>(in_c_p, data_c_p_one, in_p);
782  rst::subtract(&outi_c_p[0], &in_c_p[0], outi_c_p.size());
783  if (rst::vectorNorm(&outi_c_p[0], outi_c_p.size(), NORM_ONE) > zero) {
784  *outStream << "\n\nINCORRECT scalarMultiplyDataData (1): check scalar inverse property\n\n";
785  errorFlag = -1000;
786  }
787  art::scalarMultiplyDataData<double>(out_c_p_d, data_c_p, in_p_d);
788  art::scalarMultiplyDataData<double>(outi_c_p_d, datainv_c_p, out_c_p_d);
789  art::scalarMultiplyDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
790  rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
791  if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
792  *outStream << "\n\nINCORRECT scalarMultiplyDataData (2): check scalar inverse property\n\n";
793  errorFlag = -1000;
794  }
795  art::scalarMultiplyDataData<double>(out_c_p_d_d, data_c_p, in_p_d_d);
796  art::scalarMultiplyDataData<double>(outi_c_p_d_d, datainv_c_p, out_c_p_d_d);
797  art::scalarMultiplyDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
798  rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
799  if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
800  *outStream << "\n\nINCORRECT scalarMultiplyDataData (3): check scalar inverse property\n\n";
801  errorFlag = -1000;
802  }
803  art::scalarMultiplyDataData<double>(out_c_p_d_d, data_c_1, in_p_d_d);
804  art::scalarMultiplyDataData<double>(outi_c_p_d_d, datainv_c_1, out_c_p_d_d);
805  art::scalarMultiplyDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
806  rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
807  if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
808  *outStream << "\n\nINCORRECT scalarMultiplyDataData (4): check scalar inverse property\n\n";
809  errorFlag = -1000;
810  }
811 
812  // fill with constants
813  for (int i=0; i<in_p_d_d.size(); i++) {
814  in_p_d_d[i] = 1.0;
815  }
816  for (int i=0; i<data_c_p.size(); i++) {
817  data_c_p[i] = 5.0;
818  }
819  for (int i=0; i<data_c_1.size(); i++) {
820  data_c_1[i] = 5.0;
821  }
822  art::scalarMultiplyDataData<double>(out_c_p_d_d, data_c_p, in_p_d_d);
823  if (std::abs(rst::vectorNorm(&out_c_p_d_d[0], out_c_p_d_d.size(), NORM_ONE) - data_c_p[0]*in_p_d_d[0]*c*p*d1*d2) > zero) {
824  *outStream << "\n\nINCORRECT scalarMultiplyDataData (5): check result: "
825  << rst::vectorNorm(&out_c_p_d_d[0], out_c_p_d_d.size(), NORM_ONE) << " != "
826  << data_c_p[0]*in_p_d_d[0]*c*p*d1*d2 << "\n\n";
827  errorFlag = -1000;
828  }
829  art::scalarMultiplyDataData<double>(out_c_p_d_d, data_c_1, in_p_d_d);
830  if (std::abs(rst::vectorNorm(&out_c_p_d_d[0], out_c_p_d_d.size(), NORM_ONE) - data_c_1[0]*in_p_d_d[0]*c*p*d1*d2) > zero) {
831  *outStream << "\n\nINCORRECT scalarMultiplyDataData (6): check result: "
832  << rst::vectorNorm(&out_c_p_d_d[0], out_c_p_d_d.size(), NORM_ONE) << " != "
833  << data_c_1[0]*in_p_d_d[0]*c*p*d1*d2 << "\n\n";
834  errorFlag = -1000;
835  }
836  } // end scope
837 
838  { // start scope
839  *outStream << "\n************ Checking scalarMultiplyDataData, reciprocal=true, i.e. division (branch 1) ************\n";
840 
841  int c=5, p=9, d1=7, d2=13;
842 
843  FieldContainer<double> in_c_p(c, p);
844  FieldContainer<double> in_c_p_d(c, p, d1);
845  FieldContainer<double> in_c_p_d_d(c, p, d1, d2);
846  FieldContainer<double> data_c_p(c, p);
847  FieldContainer<double> datainv_c_p(c, p);
848  FieldContainer<double> data_c_1(c, 1);
849  FieldContainer<double> datainv_c_1(c, 1);
850  FieldContainer<double> out_c_p(c, p);
851  FieldContainer<double> outi_c_p(c, p);
852  FieldContainer<double> out_c_p_d(c, p, d1);
853  FieldContainer<double> outi_c_p_d(c, p, d1);
854  FieldContainer<double> out_c_p_d_d(c, p, d1, d2);
855  FieldContainer<double> outi_c_p_d_d(c, p, d1, d2);
856  double zero = INTREPID_TOL*10000.0;
857 
858  // fill with random numbers
859  for (int i=0; i<in_c_p.size(); i++) {
860  in_c_p[i] = Teuchos::ScalarTraits<double>::random();
861  }
862  for (int i=0; i<in_c_p_d.size(); i++) {
863  in_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
864  }
865  for (int i=0; i<in_c_p_d_d.size(); i++) {
866  in_c_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
867  }
868  for (int i=0; i<data_c_p.size(); i++) {
869  data_c_p[i] = Teuchos::ScalarTraits<double>::random();
870  datainv_c_p[i] = 1.0 / data_c_p[i];
871  }
872  for (int i=0; i<data_c_1.size(); i++) {
873  data_c_1[i] = Teuchos::ScalarTraits<double>::random();
874  datainv_c_1[i] = 1.0 / data_c_1[i];
875  }
876 
877  art::scalarMultiplyDataData<double>(out_c_p, data_c_p, in_c_p, true);
878  art::scalarMultiplyDataData<double>(outi_c_p, datainv_c_p, out_c_p, true);
879  rst::subtract(&outi_c_p[0], &in_c_p[0], outi_c_p.size());
880  if (rst::vectorNorm(&outi_c_p[0], outi_c_p.size(), NORM_ONE) > zero) {
881  *outStream << "\n\nINCORRECT scalarMultiplyDataData (1): check scalar inverse property\n\n";
882  errorFlag = -1000;
883  }
884  art::scalarMultiplyDataData<double>(out_c_p_d, data_c_p, in_c_p_d, true);
885  art::scalarMultiplyDataData<double>(outi_c_p_d, datainv_c_p, out_c_p_d, true);
886  rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
887  if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
888  *outStream << "\n\nINCORRECT scalarMultiplyDataData (2): check scalar inverse property\n\n";
889  errorFlag = -1000;
890  }
891  art::scalarMultiplyDataData<double>(out_c_p_d_d, data_c_p, in_c_p_d_d, true);
892  art::scalarMultiplyDataData<double>(outi_c_p_d_d, datainv_c_p, out_c_p_d_d, true);
893  rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
894  if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
895  *outStream << "\n\nINCORRECT scalarMultiplyDataData (3): check scalar inverse property\n\n";
896  errorFlag = -1000;
897  }
898  art::scalarMultiplyDataData<double>(out_c_p_d_d, data_c_1, in_c_p_d_d, true);
899  art::scalarMultiplyDataData<double>(outi_c_p_d_d, datainv_c_1, out_c_p_d_d, true);
900  rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
901  if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
902  *outStream << "\n\nINCORRECT scalarMultiplyDataData (4): check scalar inverse property\n\n";
903  errorFlag = -1000;
904  }
905 
906  // fill with constants
907  for (int i=0; i<in_c_p_d_d.size(); i++) {
908  in_c_p_d_d[i] = 1.0;
909  }
910  for (int i=0; i<data_c_p.size(); i++) {
911  data_c_p[i] = 5.0;
912  }
913  for (int i=0; i<data_c_1.size(); i++) {
914  data_c_1[i] = 5.0;
915  }
916  art::scalarMultiplyDataData<double>(out_c_p_d_d, data_c_p, in_c_p_d_d, true);
917  if (std::abs(rst::vectorNorm(&out_c_p_d_d[0], out_c_p_d_d.size(), NORM_ONE) - (1.0/data_c_p[0])*in_c_p_d_d[0]*c*p*d1*d2)/rst::vectorNorm(&out_c_p_d_d[0], out_c_p_d_d.size(), NORM_ONE) > zero) {
918  *outStream << "\n\nINCORRECT scalarMultiplyDataData (5): check result: "
919  << rst::vectorNorm(&out_c_p_d_d[0], out_c_p_d_d.size(), NORM_ONE) << " != "
920  << (1.0/data_c_p[0])*in_c_p_d_d[0]*c*p*d1*d2 << "\n\n";
921  errorFlag = -1000;
922  }
923  art::scalarMultiplyDataData<double>(out_c_p_d_d, data_c_1, in_c_p_d_d, true);
924  if (std::abs(rst::vectorNorm(&out_c_p_d_d[0], out_c_p_d_d.size(), NORM_ONE) - (1.0/data_c_1[0])*in_c_p_d_d[0]*c*p*d1*d2)/rst::vectorNorm(&out_c_p_d_d[0], out_c_p_d_d.size(), NORM_ONE) > zero) {
925  *outStream << "\n\nINCORRECT scalarMultiplyDataData (6): check result: "
926  << rst::vectorNorm(&out_c_p_d_d[0], out_c_p_d_d.size(), NORM_ONE) << " != "
927  << (1.0/data_c_p[0])*in_c_p_d_d[0]*c*p*d1*d2 << "\n\n";
928  errorFlag = -1000;
929  }
930  } // end scope
931 
932  { // start scope
933  *outStream << "\n************ Checking scalarMultiplyDataData, reciprocal=true, i.e. division (branch 2) ************\n";
934 
935  int c=5, p=9, d1=7, d2=13;
936 
937  FieldContainer<double> in_p(p);
938  FieldContainer<double> in_p_d(p, d1);
939  FieldContainer<double> in_p_d_d(p, d1, d2);
940  FieldContainer<double> in_c_p(c, p);
941  FieldContainer<double> in_c_p_d(c, p, d1);
942  FieldContainer<double> in_c_p_d_d(c, p, d1, d2);
943  FieldContainer<double> data_c_p(c, p);
944  FieldContainer<double> datainv_c_p(c, p);
945  FieldContainer<double> data_c_1(c, 1);
946  FieldContainer<double> datainv_c_1(c, 1);
947  FieldContainer<double> data_c_p_one(c, p);
948  FieldContainer<double> data_c_1_one(c, 1);
949  FieldContainer<double> out_c_p(c, p);
950  FieldContainer<double> outi_c_p(c, p);
951  FieldContainer<double> out_c_p_d(c, p, d1);
952  FieldContainer<double> outi_c_p_d(c, p, d1);
953  FieldContainer<double> out_c_p_d_d(c, p, d1, d2);
954  FieldContainer<double> outi_c_p_d_d(c, p, d1, d2);
955  double zero = INTREPID_TOL*10000.0;
956 
957  // fill with random numbers
958  for (int i=0; i<in_p.size(); i++) {
959  in_p[i] = Teuchos::ScalarTraits<double>::random();
960  }
961  for (int i=0; i<in_p_d.size(); i++) {
962  in_p_d[i] = Teuchos::ScalarTraits<double>::random();
963  }
964  for (int i=0; i<in_p_d_d.size(); i++) {
965  in_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
966  }
967  for (int i=0; i<data_c_p.size(); i++) {
968  data_c_p[i] = Teuchos::ScalarTraits<double>::random();
969  datainv_c_p[i] = 1.0 / data_c_p[i];
970  data_c_p_one[i] = 1.0;
971  }
972  for (int i=0; i<data_c_1.size(); i++) {
973  data_c_1[i] = Teuchos::ScalarTraits<double>::random();
974  datainv_c_1[i] = 1.0 / data_c_1[i];
975  data_c_1_one[i] = 1.0;
976  }
977 
978  art::scalarMultiplyDataData<double>(out_c_p, data_c_p, in_p, true);
979  art::scalarMultiplyDataData<double>(outi_c_p, datainv_c_p, out_c_p, true);
980  art::scalarMultiplyDataData<double>(in_c_p, data_c_p_one, in_p);
981  rst::subtract(&outi_c_p[0], &in_c_p[0], outi_c_p.size());
982  if (rst::vectorNorm(&outi_c_p[0], outi_c_p.size(), NORM_ONE) > zero) {
983  *outStream << "\n\nINCORRECT scalarMultiplyDataData (1): check scalar inverse property\n\n";
984  errorFlag = -1000;
985  }
986  art::scalarMultiplyDataData<double>(out_c_p_d, data_c_p, in_p_d, true);
987  art::scalarMultiplyDataData<double>(outi_c_p_d, datainv_c_p, out_c_p_d, true);
988  art::scalarMultiplyDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
989  rst::subtract(&outi_c_p_d[0], &in_c_p_d[0], outi_c_p_d.size());
990  if (rst::vectorNorm(&outi_c_p_d[0], outi_c_p_d.size(), NORM_ONE) > zero) {
991  *outStream << "\n\nINCORRECT scalarMultiplyDataData (2): check scalar inverse property\n\n";
992  errorFlag = -1000;
993  }
994  art::scalarMultiplyDataData<double>(out_c_p_d_d, data_c_p, in_p_d_d, true);
995  art::scalarMultiplyDataData<double>(outi_c_p_d_d, datainv_c_p, out_c_p_d_d, true);
996  art::scalarMultiplyDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
997  rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
998  if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
999  *outStream << "\n\nINCORRECT scalarMultiplyDataData (3): check scalar inverse property\n\n";
1000  errorFlag = -1000;
1001  }
1002  art::scalarMultiplyDataData<double>(out_c_p_d_d, data_c_1, in_p_d_d, true);
1003  art::scalarMultiplyDataData<double>(outi_c_p_d_d, datainv_c_1, out_c_p_d_d, true);
1004  art::scalarMultiplyDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
1005  rst::subtract(&outi_c_p_d_d[0], &in_c_p_d_d[0], outi_c_p_d_d.size());
1006  if (rst::vectorNorm(&outi_c_p_d_d[0], outi_c_p_d_d.size(), NORM_ONE) > zero) {
1007  *outStream << "\n\nINCORRECT scalarMultiplyDataData (4): check scalar inverse property\n\n";
1008  errorFlag = -1000;
1009  }
1010 
1011  // fill with constants
1012  for (int i=0; i<in_p_d_d.size(); i++) {
1013  in_p_d_d[i] = 1.0;
1014  }
1015  for (int i=0; i<data_c_p.size(); i++) {
1016  data_c_p[i] = 5.0;
1017  }
1018  for (int i=0; i<data_c_1.size(); i++) {
1019  data_c_1[i] = 5.0;
1020  }
1021  art::scalarMultiplyDataData<double>(out_c_p_d_d, data_c_p, in_p_d_d, true);
1022  if (std::abs(rst::vectorNorm(&out_c_p_d_d[0], out_c_p_d_d.size(), NORM_ONE) - (1.0/data_c_p[0])*in_p_d_d[0]*c*p*d1*d2)/rst::vectorNorm(&out_c_p_d_d[0], out_c_p_d_d.size(), NORM_ONE) > zero) {
1023  *outStream << "\n\nINCORRECT scalarMultiplyDataData (5): check result: "
1024  << rst::vectorNorm(&out_c_p_d_d[0], out_c_p_d_d.size(), NORM_ONE) << " != "
1025  << (1.0/data_c_p[0])*in_p_d_d[0]*c*p*d1*d2 << "\n\n";
1026  errorFlag = -1000;
1027  }
1028  art::scalarMultiplyDataData<double>(out_c_p_d_d, data_c_1, in_p_d_d, true);
1029  if (std::abs(rst::vectorNorm(&out_c_p_d_d[0], out_c_p_d_d.size(), NORM_ONE) - (1.0/data_c_1[0])*in_p_d_d[0]*c*p*d1*d2)/rst::vectorNorm(&out_c_p_d_d[0], out_c_p_d_d.size(), NORM_ONE) > zero) {
1030  *outStream << "\n\nINCORRECT scalarMultiplyDataData (6): check result: "
1031  << rst::vectorNorm(&out_c_p_d_d[0], out_c_p_d_d.size(), NORM_ONE) << " != "
1032  << (1.0/data_c_1[0])*in_p_d_d[0]*c*p*d1*d2 << "\n\n";
1033  errorFlag = -1000;
1034  }
1035  } // end scope
1036 
1037  /******************************************/
1038  *outStream << "\n";
1039  }
1040  catch (const std::logic_error & err) {
1041  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
1042  *outStream << err.what() << '\n';
1043  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
1044  errorFlag = -1000;
1045  };
1046 
1047 
1048  if (errorFlag != 0)
1049  std::cout << "End Result: TEST FAILED\n";
1050  else
1051  std::cout << "End Result: TEST PASSED\n";
1052 
1053  // reset format state of std::cout
1054  std::cout.copyfmt(oldFormatState);
1055 
1056  return errorFlag;
1057 }
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...