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"
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 
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 }
Intrepid::RealSpaceTools
Implementation of basic linear algebra functionality in Euclidean space.
Definition: Intrepid_RealSpaceTools.hpp:65
Intrepid_FieldContainer.hpp
Header file for utility class to provide multidimensional containers.
Intrepid_ArrayTools.hpp
Header file for utility class to provide array tools, such as tensor contractions, etc.
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:479
Intrepid_RealSpaceTools.hpp
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D...
Intrepid::ArrayTools
Utility class that provides methods for higher-order algebraic manipulation of user-defined arrays...
Definition: Intrepid_ArrayTools.hpp:67
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:55
Generated by   doxygen 1.8.5