Intrepid
test_03.cpp
Go to the documentation of this file.
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
52#include "Teuchos_oblackholestream.hpp"
53#include "Teuchos_RCP.hpp"
54#include "Teuchos_ScalarTraits.hpp"
55#include "Teuchos_GlobalMPISession.hpp"
56
57using namespace std;
58using 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
73int 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
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
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}
Header file for utility class to provide array tools, such as tensor contractions,...
Header file for utility class to provide multidimensional containers.
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D.
Utility class that provides methods for higher-order algebraic manipulation of user-defined arrays,...
static void 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...
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
Implementation of basic linear algebra functionality in Euclidean space.