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"
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 }
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 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...
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...
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...