52#include "Teuchos_oblackholestream.hpp"
53#include "Teuchos_RCP.hpp"
54#include "Teuchos_ScalarTraits.hpp"
55#include "Teuchos_GlobalMPISession.hpp"
58using namespace Intrepid;
60#define INTREPID_TEST_COMMAND( S ) \
65 catch (const std::logic_error & err) { \
66 *outStream << "Expected Error ----------------------------------------------------------------\n"; \
67 *outStream << err.what() << '\n'; \
68 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
73int main(
int argc,
char *argv[]) {
75 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
79 int iprint = argc - 1;
80 Teuchos::RCP<std::ostream> outStream;
81 Teuchos::oblackholestream bhs;
83 outStream = Teuchos::rcp(&std::cout,
false);
85 outStream = Teuchos::rcp(&bhs,
false);
88 Teuchos::oblackholestream oldFormatState;
89 oldFormatState.copyfmt(std::cout);
92 <<
"===============================================================================\n" \
94 <<
"| Unit Test (ArrayTools) |\n" \
96 <<
"| 1) Array operations: scalar multiply |\n" \
98 <<
"| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \
99 <<
"| Denis Ridzal (dridzal@sandia.gov). |\n" \
101 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
102 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
104 <<
"===============================================================================\n";
108#ifdef HAVE_INTREPID_DEBUG
109 int beginThrowNumber = Teuchos::TestForException_getThrowNumber();
110 int endThrowNumber = beginThrowNumber + 36;
114 typedef RealSpaceTools<double> rst;
115#ifdef HAVE_INTREPID_DEBUG
121 <<
"===============================================================================\n"\
122 <<
"| TEST 1: exceptions |\n"\
123 <<
"===============================================================================\n";
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);
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);
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);
149 FieldContainer<double> a_9_2(9, 2);
150 FieldContainer<double> a_10_1(10, 1);
152 FieldContainer<double> a_10_1_2(10, 1, 2);
153 FieldContainer<double> a_10_1_3(10, 1, 3);
155 FieldContainer<double> a_10_1_2_2(10, 1, 2, 2);
157 FieldContainer<double> a_2_3_2_2(2, 3, 2, 2);
158 FieldContainer<double> a_2_2_2_2(2, 2, 2, 2);
159 FieldContainer<double> a_2_10(2, 10);
160 FieldContainer<double> a_2(2);
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) );
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) );
189 FieldContainer<double> a_2_2_2(2, 2, 2);
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) );
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) );
217 catch (
const std::logic_error & err) {
218 *outStream <<
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
219 *outStream << err.what() <<
'\n';
220 *outStream <<
"-------------------------------------------------------------------------------" <<
"\n\n";
224#ifdef HAVE_INTREPID_DEBUG
225 if (Teuchos::TestForException_getThrowNumber() != endThrowNumber)
231 <<
"===============================================================================\n"\
232 <<
"| TEST 2: correctness of math operations |\n"\
233 <<
"===============================================================================\n";
235 outStream->precision(20);
239 *outStream <<
"\n************ Checking scalarMultiplyDataField, reciprocal=false (branch 1) ************\n";
241 int c=5, p=9, f=7, d1=7, d2=13;
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;
259 for (
int i=0; i<in_c_f_p.size(); i++) {
260 in_c_f_p[i] = Teuchos::ScalarTraits<double>::random();
262 for (
int i=0; i<in_c_f_p_d.size(); i++) {
263 in_c_f_p_d[i] = Teuchos::ScalarTraits<double>::random();
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();
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];
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];
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";
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";
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";
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";
307 for (
int i=0; i<in_c_f_p_d_d.size(); i++) {
308 in_c_f_p_d_d[i] = 1.0;
310 for (
int i=0; i<data_c_p.size(); i++) {
313 for (
int i=0; i<data_c_1.size(); i++) {
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";
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";
333 *outStream <<
"\n************ Checking scalarMultiplyDataField, reciprocal=false (branch 2) ************\n";
335 int c=5, p=9, f=7, d1=7, d2=13;
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;
358 for (
int i=0; i<in_f_p.size(); i++) {
359 in_f_p[i] = Teuchos::ScalarTraits<double>::random();
361 for (
int i=0; i<in_f_p_d.size(); i++) {
362 in_f_p_d[i] = Teuchos::ScalarTraits<double>::random();
364 for (
int i=0; i<in_f_p_d_d.size(); i++) {
365 in_f_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
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;
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;
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";
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";
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";
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";
412 for (
int i=0; i<in_f_p_d_d.size(); i++) {
415 for (
int i=0; i<data_c_p.size(); i++) {
418 for (
int i=0; i<data_c_1.size(); i++) {
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";
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";
438 *outStream <<
"\n************ Checking scalarMultiplyDataField, reciprocal=true, i.e. division (branch 1) ************\n";
440 int c=5, p=9, f=7, d1=7, d2=13;
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;
458 for (
int i=0; i<in_c_f_p.size(); i++) {
459 in_c_f_p[i] = Teuchos::ScalarTraits<double>::random();
461 for (
int i=0; i<in_c_f_p_d.size(); i++) {
462 in_c_f_p_d[i] = Teuchos::ScalarTraits<double>::random();
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();
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];
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];
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";
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";
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";
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";
506 for (
int i=0; i<in_c_f_p_d_d.size(); i++) {
507 in_c_f_p_d_d[i] = 1.0;
509 for (
int i=0; i<data_c_p.size(); i++) {
512 for (
int i=0; i<data_c_1.size(); i++) {
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";
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";
532 *outStream <<
"\n************ Checking scalarMultiplyDataField, reciprocal=true, i.e. division (branch 2) ************\n";
534 int c=5, p=9, f=7, d1=7, d2=13;
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;
557 for (
int i=0; i<in_f_p.size(); i++) {
558 in_f_p[i] = Teuchos::ScalarTraits<double>::random();
560 for (
int i=0; i<in_f_p_d.size(); i++) {
561 in_f_p_d[i] = Teuchos::ScalarTraits<double>::random();
563 for (
int i=0; i<in_f_p_d_d.size(); i++) {
564 in_f_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
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;
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;
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";
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";
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";
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";
611 for (
int i=0; i<in_f_p_d_d.size(); i++) {
614 for (
int i=0; i<data_c_p.size(); i++) {
617 for (
int i=0; i<data_c_1.size(); i++) {
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";
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";
640 *outStream <<
"\n************ Checking scalarMultiplyDataData, reciprocal=false (branch 1) ************\n";
642 int c=5, p=9, d1=7, d2=13;
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;
660 for (
int i=0; i<in_c_p.size(); i++) {
661 in_c_p[i] = Teuchos::ScalarTraits<double>::random();
663 for (
int i=0; i<in_c_p_d.size(); i++) {
664 in_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
666 for (
int i=0; i<in_c_p_d_d.size(); i++) {
667 in_c_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
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];
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];
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";
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";
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";
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";
708 for (
int i=0; i<in_c_p_d_d.size(); i++) {
711 for (
int i=0; i<data_c_p.size(); i++) {
714 for (
int i=0; i<data_c_1.size(); i++) {
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";
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";
734 *outStream <<
"\n************ Checking scalarMultiplyDataData, reciprocal=false (branch 2) ************\n";
736 int c=5, p=9, d1=7, d2=13;
738 FieldContainer<double> in_p(p);
739 FieldContainer<double> in_p_d(p, d1);
740 FieldContainer<double> in_p_d_d(p, d1, d2);
741 FieldContainer<double> in_c_p(c, p);
742 FieldContainer<double> in_c_p_d(c, p, d1);
743 FieldContainer<double> in_c_p_d_d(c, p, d1, d2);
744 FieldContainer<double> data_c_p(c, p);
745 FieldContainer<double> datainv_c_p(c, p);
746 FieldContainer<double> data_c_1(c, 1);
747 FieldContainer<double> datainv_c_1(c, 1);
748 FieldContainer<double> data_c_p_one(c, p);
749 FieldContainer<double> data_c_1_one(c, 1);
750 FieldContainer<double> out_c_p(c, p);
751 FieldContainer<double> outi_c_p(c, p);
752 FieldContainer<double> out_c_p_d(c, p, d1);
753 FieldContainer<double> outi_c_p_d(c, p, d1);
754 FieldContainer<double> out_c_p_d_d(c, p, d1, d2);
755 FieldContainer<double> outi_c_p_d_d(c, p, d1, d2);
756 double zero = INTREPID_TOL*10000.0;
759 for (
int i=0; i<in_p.size(); i++) {
760 in_p[i] = Teuchos::ScalarTraits<double>::random();
762 for (
int i=0; i<in_p_d.size(); i++) {
763 in_p_d[i] = Teuchos::ScalarTraits<double>::random();
765 for (
int i=0; i<in_p_d_d.size(); i++) {
766 in_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
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;
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;
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";
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";
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";
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";
813 for (
int i=0; i<in_p_d_d.size(); i++) {
816 for (
int i=0; i<data_c_p.size(); i++) {
819 for (
int i=0; i<data_c_1.size(); i++) {
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";
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";
839 *outStream <<
"\n************ Checking scalarMultiplyDataData, reciprocal=true, i.e. division (branch 1) ************\n";
841 int c=5, p=9, d1=7, d2=13;
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;
859 for (
int i=0; i<in_c_p.size(); i++) {
860 in_c_p[i] = Teuchos::ScalarTraits<double>::random();
862 for (
int i=0; i<in_c_p_d.size(); i++) {
863 in_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
865 for (
int i=0; i<in_c_p_d_d.size(); i++) {
866 in_c_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
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];
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];
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";
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";
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";
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";
907 for (
int i=0; i<in_c_p_d_d.size(); i++) {
910 for (
int i=0; i<data_c_p.size(); i++) {
913 for (
int i=0; i<data_c_1.size(); i++) {
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";
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";
933 *outStream <<
"\n************ Checking scalarMultiplyDataData, reciprocal=true, i.e. division (branch 2) ************\n";
935 int c=5, p=9, d1=7, d2=13;
937 FieldContainer<double> in_p(p);
938 FieldContainer<double> in_p_d(p, d1);
939 FieldContainer<double> in_p_d_d(p, d1, d2);
940 FieldContainer<double> in_c_p(c, p);
941 FieldContainer<double> in_c_p_d(c, p, d1);
942 FieldContainer<double> in_c_p_d_d(c, p, d1, d2);
943 FieldContainer<double> data_c_p(c, p);
944 FieldContainer<double> datainv_c_p(c, p);
945 FieldContainer<double> data_c_1(c, 1);
946 FieldContainer<double> datainv_c_1(c, 1);
947 FieldContainer<double> data_c_p_one(c, p);
948 FieldContainer<double> data_c_1_one(c, 1);
949 FieldContainer<double> out_c_p(c, p);
950 FieldContainer<double> outi_c_p(c, p);
951 FieldContainer<double> out_c_p_d(c, p, d1);
952 FieldContainer<double> outi_c_p_d(c, p, d1);
953 FieldContainer<double> out_c_p_d_d(c, p, d1, d2);
954 FieldContainer<double> outi_c_p_d_d(c, p, d1, d2);
955 double zero = INTREPID_TOL*10000.0;
958 for (
int i=0; i<in_p.size(); i++) {
959 in_p[i] = Teuchos::ScalarTraits<double>::random();
961 for (
int i=0; i<in_p_d.size(); i++) {
962 in_p_d[i] = Teuchos::ScalarTraits<double>::random();
964 for (
int i=0; i<in_p_d_d.size(); i++) {
965 in_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
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;
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;
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";
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";
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";
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";
1012 for (
int i=0; i<in_p_d_d.size(); i++) {
1015 for (
int i=0; i<data_c_p.size(); i++) {
1018 for (
int i=0; i<data_c_1.size(); i++) {
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";
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";
1040 catch (
const std::logic_error & err) {
1041 *outStream <<
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
1042 *outStream << err.what() <<
'\n';
1043 *outStream <<
"-------------------------------------------------------------------------------" <<
"\n\n";
1049 std::cout <<
"End Result: TEST FAILED\n";
1051 std::cout <<
"End Result: TEST PASSED\n";
1054 std::cout.copyfmt(oldFormatState);
Header file for utility class to provide multidimensional containers.