Intrepid
test_03.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: dot 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);
129  FieldContainer<double> a_2_2(2, 2);
130  FieldContainer<double> a_3_2(3, 2);
131  FieldContainer<double> a_2_2_2(2, 2, 2);
132  FieldContainer<double> a_2_2_2_2(2, 2, 2, 2);
133  FieldContainer<double> a_10_1(10, 1);
134  FieldContainer<double> a_10_2(10, 2);
135  FieldContainer<double> a_10_3(10, 3);
136  FieldContainer<double> a_10_1_2(10, 1, 2);
137  FieldContainer<double> a_10_2_2(10, 2, 2);
138  FieldContainer<double> a_10_2_2_2(10, 2, 2, 2);
139  FieldContainer<double> a_10_2_2_2_2(10, 2, 2, 2, 2);
140  FieldContainer<double> a_9_2_2(9, 2, 2);
141  FieldContainer<double> a_10_3_2(10, 3, 2);
142  FieldContainer<double> a_10_2_3(10, 2, 3);
143  FieldContainer<double> a_10_2_2_3(10, 2, 2, 3);
144  FieldContainer<double> a_10_1_2_2(10, 1, 2, 2);
145 
146  *outStream << "-> dotMultiplyDataField:\n";
147  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_2_2, a_2, a_2_2) );
148  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_2_2, a_10_2_2, a_10_2_2_2_2) );
149  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2, a_10_2_2, a_10_2_2_2) );
150  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_9_2_2, a_10_2_2_2) );
151  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_10_3_2, a_10_2_2_2) );
152  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_10_2_3, a_10_2_2_2) );
153  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_10_2_2_3, a_10_2_2_2_2) );
154  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_2_2_2, a_10_2_2_2, a_10_2_2_2_2) );
155  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_3_2, a_10_2_2_2, a_10_2_2_2_2) );
156  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_3, a_10_2_2_2, a_10_2_2_2_2) );
157  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_10_2_2_2, a_10_2_2_2_2) );
158  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_10_1_2_2, a_10_2_2_2_2) );
159  //
160  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_2_2, a_2, a_2) );
161  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_2_2, a_10_2_2, a_10_2) );
162  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2, a_10_2_2, a_10_2_2) );
163  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_10_2_2, a_10_3_2) );
164  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_10_2_2, a_10_2_2) );
165  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_3, a_10_2_2, a_2_2_2) );
166  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_9_2_2, a_2_2_2) );
167  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_10_2_3, a_2_2_2) );
168  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_10_2_2_3, a_2_2_2_2) );
169  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_10_2_2_2, a_2_2_2_2) );
170  INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_10_1, a_2_2) );
171 
172  *outStream << "-> dotMultiplyDataData:\n";
173  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_2_2, a_2, a_2_2) );
174  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_2_2, a_10_2_2, a_10_2_2_2) );
175  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2_2, a_10_2_2, a_10_2_2) );
176  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_9_2_2, a_10_2_2) );
177  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_3_2, a_10_2_2) );
178  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_2_3, a_10_2_2) );
179  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_2_2_2, a_10_2_2_3) );
180  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_9_2_2, a_9_2_2) );
181  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_3_2, a_10_3_2) );
182  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_2, a_10_2) );
183  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_2_2, a_10_2_2) );
184  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_2_2_2, a_10_2_2_2) );
185  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_1, a_10_2) );
186  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_1_2, a_10_2_2) );
187  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_1_2_2, a_10_2_2_2) );
188  //
189  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_2_2, a_10_2_2_2_2, a_10_2) );
190  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_2_2, a_10_2_2, a_2) );
191  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2_2, a_10_2_2, a_10_2) );
192  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_2_2, a_10_2) );
193  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_3, a_10_2_2, a_2_2) );
194  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_9_2_2, a_2_2) );
195  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_2_3, a_2_2) );
196  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_2_2_3, a_2_2_2) );
197  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_2, a_2) );
198  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_2_2, a_2_2) );
199  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_2_2_2, a_2_2_2) );
200  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_1, a_2) );
201  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_1_2, a_2_2) );
202  INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_1_2_2, a_2_2_2) );
203 #endif
204 
205  }
206  catch (const std::logic_error & err) {
207  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
208  *outStream << err.what() << '\n';
209  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
210  errorFlag = -1000;
211  };
212 
213 #ifdef HAVE_INTREPID_DEBUG
214  if (Teuchos::TestForException_getThrowNumber() != endThrowNumber)
215  errorFlag++;
216 #endif
217 
218  *outStream \
219  << "\n"
220  << "===============================================================================\n"\
221  << "| TEST 2: correctness of math operations |\n"\
222  << "===============================================================================\n";
223 
224  outStream->precision(20);
225 
226  try {
227  { // start scope
228  *outStream << "\n************ Checking dotMultiplyDataField, (branch 1) ************\n";
229 
230  int c=5, p=9, f=7, d1=6, d2=14;
231 
232  FieldContainer<double> in_c_f_p(c, f, p);
233  FieldContainer<double> in_c_f_p_d(c, f, p, d1);
234  FieldContainer<double> in_c_f_p_d_d(c, f, p, d1, d2);
235  FieldContainer<double> data_c_p(c, p);
236  FieldContainer<double> data_c_p_d(c, p, d1);
237  FieldContainer<double> data_c_p_d_d(c, p, d1, d2);
238  FieldContainer<double> data_c_1(c, 1);
239  FieldContainer<double> data_c_1_d(c, 1, d1);
240  FieldContainer<double> data_c_1_d_d(c, 1, d1, d2);
241  FieldContainer<double> out_c_f_p(c, f, p);
242  FieldContainer<double> outSM_c_f_p(c, f, p);
243  FieldContainer<double> outDM_c_f_p(c, f, p);
244  double zero = INTREPID_TOL*10000.0;
245 
246  // fill with random numbers
247  for (int i=0; i<in_c_f_p.size(); i++) {
248  in_c_f_p[i] = Teuchos::ScalarTraits<double>::random();
249  }
250  // fill with alternating 1's and -1's
251  for (int i=0; i<in_c_f_p_d.size(); i++) {
252  in_c_f_p_d[i] = std::pow((double)-1.0, (int)i);
253  }
254  // fill with 1's
255  for (int i=0; i<in_c_f_p_d_d.size(); i++) {
256  in_c_f_p_d_d[i] = 1.0;
257  }
258  // fill with random numbers
259  for (int i=0; i<data_c_p.size(); i++) {
260  data_c_p[i] = Teuchos::ScalarTraits<double>::random();
261  }
262  // fill with 1's
263  for (int i=0; i<data_c_1.size(); i++) {
264  data_c_1[i] = Teuchos::ScalarTraits<double>::random();
265  }
266  // fill with 1's
267  for (int i=0; i<data_c_p_d.size(); i++) {
268  data_c_p_d[i] = 1.0;
269  }
270  // fill with 1's
271  for (int i=0; i<data_c_1_d.size(); i++) {
272  data_c_1_d[i] = 1.0;
273  }
274  // fill with 1's
275  for (int i=0; i<data_c_p_d_d.size(); i++) {
276  data_c_p_d_d[i] = 1.0;
277  }
278  // fill with 1's
279  for (int i=0; i<data_c_1_d_d.size(); i++) {
280  data_c_1_d_d[i] = 1.0;
281  }
282 
283  art::scalarMultiplyDataField<double>(outSM_c_f_p, data_c_p, in_c_f_p);
284  art::dotMultiplyDataField<double>(outDM_c_f_p, data_c_p, in_c_f_p);
285  rst::subtract(out_c_f_p, outSM_c_f_p, outDM_c_f_p);
286  if (rst::vectorNorm(&out_c_f_p[0], out_c_f_p.size(), NORM_ONE) > zero) {
287  *outStream << "\n\nINCORRECT dotMultiplyDataField (1): check dot multiply for scalars vs. scalar multiply\n\n";
288  errorFlag = -1000;
289  }
290  art::dotMultiplyDataField<double>(out_c_f_p, data_c_p_d, in_c_f_p_d);
291  if (rst::vectorNorm(&out_c_f_p[0], out_c_f_p.size(), NORM_ONE) > zero) {
292  *outStream << "\n\nINCORRECT dotMultiplyDataField (2): check dot multiply of orthogonal vectors\n\n";
293  errorFlag = -1000;
294  }
295  art::dotMultiplyDataField<double>(out_c_f_p, data_c_p_d_d, in_c_f_p_d_d);
296  if ((rst::vectorNorm(&out_c_f_p[0], out_c_f_p.size(), NORM_INF) - d1*d2) > zero) {
297  *outStream << "\n\nINCORRECT dotMultiplyDataField (3): check dot multiply for tensors of 1s\n\n";
298  errorFlag = -1000;
299  }
300  art::scalarMultiplyDataField<double>(outSM_c_f_p, data_c_1, in_c_f_p);
301  art::dotMultiplyDataField<double>(outDM_c_f_p, data_c_1, in_c_f_p);
302  rst::subtract(out_c_f_p, outSM_c_f_p, outDM_c_f_p);
303  if (rst::vectorNorm(&out_c_f_p[0], out_c_f_p.size(), NORM_ONE) > zero) {
304  *outStream << "\n\nINCORRECT dotMultiplyDataField (4): check dot multiply for scalars vs. scalar multiply\n\n";
305  errorFlag = -1000;
306  }
307  art::dotMultiplyDataField<double>(out_c_f_p, data_c_1_d, in_c_f_p_d);
308  if (rst::vectorNorm(&out_c_f_p[0], out_c_f_p.size(), NORM_ONE) > zero) {
309  *outStream << "\n\nINCORRECT dotMultiplyDataField (5): check dot multiply of orthogonal vectors\n\n";
310  errorFlag = -1000;
311  }
312  art::dotMultiplyDataField<double>(out_c_f_p, data_c_1_d_d, in_c_f_p_d_d);
313  if ((rst::vectorNorm(&out_c_f_p[0], out_c_f_p.size(), NORM_INF) - d1*d2) > zero) {
314  *outStream << "\n\nINCORRECT dotMultiplyDataField (6): check dot multiply for tensors of 1s\n\n";
315  errorFlag = -1000;
316  }
317  } // end scope
318 
319 
320  { // start scope
321  *outStream << "\n************ Checking dotMultiplyDataField, (branch 2) ************\n";
322 
323  int c=5, p=9, f=7, d1=6, d2=14;
324 
325  FieldContainer<double> in_f_p(f, p);
326  FieldContainer<double> in_f_p_d(f, p, d1);
327  FieldContainer<double> in_f_p_d_d(f, p, d1, d2);
328  FieldContainer<double> data_c_p(c, p);
329  FieldContainer<double> data_c_p_d(c, p, d1);
330  FieldContainer<double> data_c_p_d_d(c, p, d1, d2);
331  FieldContainer<double> data_c_1(c, 1);
332  FieldContainer<double> data_c_1_d(c, 1, d1);
333  FieldContainer<double> data_c_1_d_d(c, 1, d1, d2);
334  FieldContainer<double> out_c_f_p(c, f, p);
335  FieldContainer<double> outSM_c_f_p(c, f, p);
336  FieldContainer<double> outDM_c_f_p(c, f, p);
337  double zero = INTREPID_TOL*10000.0;
338 
339  // fill with random numbers
340  for (int i=0; i<in_f_p.size(); i++) {
341  in_f_p[i] = Teuchos::ScalarTraits<double>::random();
342  }
343  // fill with alternating 1's and -1's
344  for (int i=0; i<in_f_p_d.size(); i++) {
345  in_f_p_d[i] = std::pow((double)-1.0, (int)i);
346  }
347  // fill with 1's
348  for (int i=0; i<in_f_p_d_d.size(); i++) {
349  in_f_p_d_d[i] = 1.0;
350  }
351  // fill with random numbers
352  for (int i=0; i<data_c_p.size(); i++) {
353  data_c_p[i] = Teuchos::ScalarTraits<double>::random();
354  }
355  // fill with 1's
356  for (int i=0; i<data_c_1.size(); i++) {
357  data_c_1[i] = Teuchos::ScalarTraits<double>::random();
358  }
359  // fill with 1's
360  for (int i=0; i<data_c_p_d.size(); i++) {
361  data_c_p_d[i] = 1.0;
362  }
363  // fill with 1's
364  for (int i=0; i<data_c_1_d.size(); i++) {
365  data_c_1_d[i] = 1.0;
366  }
367  // fill with 1's
368  for (int i=0; i<data_c_p_d_d.size(); i++) {
369  data_c_p_d_d[i] = 1.0;
370  }
371  // fill with 1's
372  for (int i=0; i<data_c_1_d_d.size(); i++) {
373  data_c_1_d_d[i] = 1.0;
374  }
375 
376  art::scalarMultiplyDataField<double>(outSM_c_f_p, data_c_p, in_f_p);
377  art::dotMultiplyDataField<double>(outDM_c_f_p, data_c_p, in_f_p);
378  rst::subtract(out_c_f_p, outSM_c_f_p, outDM_c_f_p);
379  if (rst::vectorNorm(&out_c_f_p[0], out_c_f_p.size(), NORM_ONE) > zero) {
380  *outStream << "\n\nINCORRECT dotMultiplyDataField (7): check dot multiply for scalars vs. scalar multiply\n\n";
381  errorFlag = -1000;
382  }
383  art::dotMultiplyDataField<double>(out_c_f_p, data_c_p_d, in_f_p_d);
384  if (rst::vectorNorm(&out_c_f_p[0], out_c_f_p.size(), NORM_ONE) > zero) {
385  *outStream << "\n\nINCORRECT dotMultiplyDataField (8): check dot multiply of orthogonal vectors\n\n";
386  errorFlag = -1000;
387  }
388  art::dotMultiplyDataField<double>(out_c_f_p, data_c_p_d_d, in_f_p_d_d);
389  if ((rst::vectorNorm(&out_c_f_p[0], out_c_f_p.size(), NORM_INF) - d1*d2) > zero) {
390  *outStream << "\n\nINCORRECT dotMultiplyDataField (9): check dot multiply for tensors of 1s\n\n";
391  errorFlag = -1000;
392  }
393  art::scalarMultiplyDataField<double>(outSM_c_f_p, data_c_1, in_f_p);
394  art::dotMultiplyDataField<double>(outDM_c_f_p, data_c_1, in_f_p);
395  rst::subtract(out_c_f_p, outSM_c_f_p, outDM_c_f_p);
396  if (rst::vectorNorm(&out_c_f_p[0], out_c_f_p.size(), NORM_ONE) > zero) {
397  *outStream << "\n\nINCORRECT dotMultiplyDataField (10): check dot multiply for scalars vs. scalar multiply\n\n";
398  errorFlag = -1000;
399  }
400  art::dotMultiplyDataField<double>(out_c_f_p, data_c_1_d, in_f_p_d);
401  if (rst::vectorNorm(&out_c_f_p[0], out_c_f_p.size(), NORM_ONE) > zero) {
402  *outStream << "\n\nINCORRECT dotMultiplyDataField (11): check dot multiply of orthogonal vectors\n\n";
403  errorFlag = -1000;
404  }
405  art::dotMultiplyDataField<double>(out_c_f_p, data_c_1_d_d, in_f_p_d_d);
406  if ((rst::vectorNorm(&out_c_f_p[0], out_c_f_p.size(), NORM_INF) - d1*d2) > zero) {
407  *outStream << "\n\nINCORRECT dotMultiplyDataField (12): check dot multiply for tensors of 1s\n\n";
408  errorFlag = -1000;
409  }
410  } // end scope
411 
412 
413 
414 
415  { // start scope
416  *outStream << "\n************ Checking dotMultiplyDataData, (branch 1) ************\n";
417 
418  int c=5, p=9, d1=6, d2=14;
419 
420  FieldContainer<double> in_c_p(c, p);
421  FieldContainer<double> in_c_p_d(c, p, d1);
422  FieldContainer<double> in_c_p_d_d(c, p, d1, d2);
423  FieldContainer<double> data_c_p(c, p);
424  FieldContainer<double> data_c_p_d(c, p, d1);
425  FieldContainer<double> data_c_p_d_d(c, p, d1, d2);
426  FieldContainer<double> data_c_1(c, 1);
427  FieldContainer<double> data_c_1_d(c, 1, d1);
428  FieldContainer<double> data_c_1_d_d(c, 1, d1, d2);
429  FieldContainer<double> out_c_p(c, p);
430  FieldContainer<double> outSM_c_p(c, p);
431  FieldContainer<double> outDM_c_p(c, p);
432  double zero = INTREPID_TOL*10000.0;
433 
434  // fill with random numbers
435  for (int i=0; i<in_c_p.size(); i++) {
436  in_c_p[i] = Teuchos::ScalarTraits<double>::random();
437  }
438  // fill with alternating 1's and -1's
439  for (int i=0; i<in_c_p_d.size(); i++) {
440  in_c_p_d[i] = std::pow((double)-1.0, (int)i);
441  }
442  // fill with 1's
443  for (int i=0; i<in_c_p_d_d.size(); i++) {
444  in_c_p_d_d[i] = 1.0;
445  }
446  // fill with random numbers
447  for (int i=0; i<data_c_p.size(); i++) {
448  data_c_p[i] = Teuchos::ScalarTraits<double>::random();
449  }
450  // fill with 1's
451  for (int i=0; i<data_c_1.size(); i++) {
452  data_c_1[i] = Teuchos::ScalarTraits<double>::random();
453  }
454  // fill with 1's
455  for (int i=0; i<data_c_p_d.size(); i++) {
456  data_c_p_d[i] = 1.0;
457  }
458  // fill with 1's
459  for (int i=0; i<data_c_1_d.size(); i++) {
460  data_c_1_d[i] = 1.0;
461  }
462  // fill with 1's
463  for (int i=0; i<data_c_p_d_d.size(); i++) {
464  data_c_p_d_d[i] = 1.0;
465  }
466  // fill with 1's
467  for (int i=0; i<data_c_1_d_d.size(); i++) {
468  data_c_1_d_d[i] = 1.0;
469  }
470 
471  art::scalarMultiplyDataData<double>(outSM_c_p, data_c_p, in_c_p);
472  art::dotMultiplyDataData<double>(outDM_c_p, data_c_p, in_c_p);
473  rst::subtract(out_c_p, outSM_c_p, outDM_c_p);
474  if (rst::vectorNorm(&out_c_p[0], out_c_p.size(), NORM_ONE) > zero) {
475  *outStream << "\n\nINCORRECT dotMultiplyDataData (1): check dot multiply for scalars vs. scalar multiply\n\n";
476  errorFlag = -1000;
477  }
478  art::dotMultiplyDataData<double>(out_c_p, data_c_p_d, in_c_p_d);
479  if (rst::vectorNorm(&out_c_p[0], out_c_p.size(), NORM_ONE) > zero) {
480  *outStream << "\n\nINCORRECT dotMultiplyDataData (2): check dot multiply of orthogonal vectors\n\n";
481  errorFlag = -1000;
482  }
483  art::dotMultiplyDataData<double>(out_c_p, data_c_p_d_d, in_c_p_d_d);
484  if ((rst::vectorNorm(&out_c_p[0], out_c_p.size(), NORM_INF) - d1*d2) > zero) {
485  *outStream << "\n\nINCORRECT dotMultiplyDataData (3): check dot multiply for tensors of 1s\n\n";
486  errorFlag = -1000;
487  }
488  art::scalarMultiplyDataData<double>(outSM_c_p, data_c_1, in_c_p);
489  art::dotMultiplyDataData<double>(outDM_c_p, data_c_1, in_c_p);
490  rst::subtract(out_c_p, outSM_c_p, outDM_c_p);
491  if (rst::vectorNorm(&out_c_p[0], out_c_p.size(), NORM_ONE) > zero) {
492  *outStream << "\n\nINCORRECT dotMultiplyDataData (4): check dot multiply for scalars vs. scalar multiply\n\n";
493  errorFlag = -1000;
494  }
495  art::dotMultiplyDataData<double>(out_c_p, data_c_1_d, in_c_p_d);
496  if (rst::vectorNorm(&out_c_p[0], out_c_p.size(), NORM_ONE) > zero) {
497  *outStream << "\n\nINCORRECT dotMultiplyDataData (5): check dot multiply of orthogonal vectors\n\n";
498  errorFlag = -1000;
499  }
500  art::dotMultiplyDataData<double>(out_c_p, data_c_1_d_d, in_c_p_d_d);
501  if ((rst::vectorNorm(&out_c_p[0], out_c_p.size(), NORM_INF) - d1*d2) > zero) {
502  *outStream << "\n\nINCORRECT dotMultiplyDataData (6): check dot multiply for tensors of 1s\n\n";
503  errorFlag = -1000;
504  }
505  } // end scope
506 
507 
508  { // start scope
509  *outStream << "\n************ Checking dotMultiplyDataData, (branch 2) ************\n";
510 
511  int c=5, p=9, d1=6, d2=14;
512 
513  FieldContainer<double> in_p(p);
514  FieldContainer<double> in_p_d(p, d1);
515  FieldContainer<double> in_p_d_d(p, d1, d2);
516  FieldContainer<double> data_c_p(c, p);
517  FieldContainer<double> data_c_p_d(c, p, d1);
518  FieldContainer<double> data_c_p_d_d(c, p, d1, d2);
519  FieldContainer<double> data_c_1(c, 1);
520  FieldContainer<double> data_c_1_d(c, 1, d1);
521  FieldContainer<double> data_c_1_d_d(c, 1, d1, d2);
522  FieldContainer<double> out_c_p(c, p);
523  FieldContainer<double> outSM_c_p(c, p);
524  FieldContainer<double> outDM_c_p(c, p);
525  double zero = INTREPID_TOL*10000.0;
526 
527  // fill with random numbers
528  for (int i=0; i<in_p.size(); i++) {
529  in_p[i] = Teuchos::ScalarTraits<double>::random();
530  }
531  // fill with alternating 1's and -1's
532  for (int i=0; i<in_p_d.size(); i++) {
533  in_p_d[i] = std::pow((double)-1.0, (int)i);
534  }
535  // fill with 1's
536  for (int i=0; i<in_p_d_d.size(); i++) {
537  in_p_d_d[i] = 1.0;
538  }
539  // fill with random numbers
540  for (int i=0; i<data_c_p.size(); i++) {
541  data_c_p[i] = Teuchos::ScalarTraits<double>::random();
542  }
543  // fill with 1's
544  for (int i=0; i<data_c_1.size(); i++) {
545  data_c_1[i] = Teuchos::ScalarTraits<double>::random();
546  }
547  // fill with 1's
548  for (int i=0; i<data_c_p_d.size(); i++) {
549  data_c_p_d[i] = 1.0;
550  }
551  // fill with 1's
552  for (int i=0; i<data_c_1_d.size(); i++) {
553  data_c_1_d[i] = 1.0;
554  }
555  // fill with 1's
556  for (int i=0; i<data_c_p_d_d.size(); i++) {
557  data_c_p_d_d[i] = 1.0;
558  }
559  // fill with 1's
560  for (int i=0; i<data_c_1_d_d.size(); i++) {
561  data_c_1_d_d[i] = 1.0;
562  }
563 
564  art::scalarMultiplyDataData<double>(outSM_c_p, data_c_p, in_p);
565  art::dotMultiplyDataData<double>(outDM_c_p, data_c_p, in_p);
566  rst::subtract(out_c_p, outSM_c_p, outDM_c_p);
567  if (rst::vectorNorm(&out_c_p[0], out_c_p.size(), NORM_ONE) > zero) {
568  *outStream << "\n\nINCORRECT dotMultiplyDataData (7): check dot multiply for scalars vs. scalar multiply\n\n";
569  errorFlag = -1000;
570  }
571  art::dotMultiplyDataData<double>(out_c_p, data_c_p_d, in_p_d);
572  if (rst::vectorNorm(&out_c_p[0], out_c_p.size(), NORM_ONE) > zero) {
573  *outStream << "\n\nINCORRECT dotMultiplyDataData (8): check dot multiply of orthogonal vectors\n\n";
574  errorFlag = -1000;
575  }
576  art::dotMultiplyDataData<double>(out_c_p, data_c_p_d_d, in_p_d_d);
577  if ((rst::vectorNorm(&out_c_p[0], out_c_p.size(), NORM_INF) - d1*d2) > zero) {
578  *outStream << "\n\nINCORRECT dotMultiplyDataData (9): check dot multiply for tensors of 1s\n\n";
579  errorFlag = -1000;
580  }
581  art::scalarMultiplyDataData<double>(outSM_c_p, data_c_1, in_p);
582  art::dotMultiplyDataData<double>(outDM_c_p, data_c_1, in_p);
583  rst::subtract(out_c_p, outSM_c_p, outDM_c_p);
584  if (rst::vectorNorm(&out_c_p[0], out_c_p.size(), NORM_ONE) > zero) {
585  *outStream << "\n\nINCORRECT dotMultiplyDataData (10): check dot multiply for scalars vs. scalar multiply\n\n";
586  errorFlag = -1000;
587  }
588  art::dotMultiplyDataData<double>(out_c_p, data_c_1_d, in_p_d);
589  if (rst::vectorNorm(&out_c_p[0], out_c_p.size(), NORM_ONE) > zero) {
590  *outStream << "\n\nINCORRECT dotMultiplyDataData (11): check dot multiply of orthogonal vectors\n\n";
591  errorFlag = -1000;
592  }
593  art::dotMultiplyDataData<double>(out_c_p, data_c_1_d_d, in_p_d_d);
594  if ((rst::vectorNorm(&out_c_p[0], out_c_p.size(), NORM_INF) - d1*d2) > zero) {
595  *outStream << "\n\nINCORRECT dotMultiplyDataData (12): check dot multiply for tensors of 1s\n\n";
596  errorFlag = -1000;
597  }
598  } // end scope
599 
600  /******************************************/
601  *outStream << "\n";
602  }
603  catch (const std::logic_error & err) {
604  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
605  *outStream << err.what() << '\n';
606  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
607  errorFlag = -1000;
608  };
609 
610 
611  if (errorFlag != 0)
612  std::cout << "End Result: TEST FAILED\n";
613  else
614  std::cout << "End Result: TEST PASSED\n";
615 
616  // reset format state of std::cout
617  std::cout.copyfmt(oldFormatState);
618 
619  return errorFlag;
620 }
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::dotMultiplyDataField
static void dotMultiplyDataField(ArrayOutFields &outputFields, const ArrayInData &inputDataLeft, const ArrayInFields &inputFields)
There are two use cases: (1) dot product of a rank-3, 4 or 5 container inputFields with dimensions (C...
Definition: Intrepid_ArrayToolsDefDot.hpp:52
Intrepid::ArrayTools::dotMultiplyDataData
static void dotMultiplyDataData(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight)
There are two use cases: (1) dot product of a rank-2, 3 or 4 container inputDataRight with dimensions...
Definition: Intrepid_ArrayToolsDefDot.hpp:351
Intrepid_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
Generated by   doxygen 1.8.5