Intrepid
test_02.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: 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);
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
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
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}
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 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...
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,...
Implementation of a templated lexicographical container for a multi-indexed scalar quantity....
Implementation of basic linear algebra functionality in Euclidean space.