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: contractions |\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 + 81;
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 *outStream <<
"-> contractFieldFieldScalar:\n";
137 INTREPID_TEST_COMMAND( atools.contractFieldFieldScalar<
double>(a_2_2, a_2_2, a_2_2, COMP_CPP) );
138 INTREPID_TEST_COMMAND( atools.contractFieldFieldScalar<
double>(a_2_2, a_10_2_2, a_2_2, COMP_CPP) );
139 INTREPID_TEST_COMMAND( atools.contractFieldFieldScalar<
double>(a_2_2, a_10_2_2, a_10_2_2, COMP_CPP) );
140 INTREPID_TEST_COMMAND( atools.contractFieldFieldScalar<
double>(a_10_2_2, a_9_2_2, a_10_2_2, COMP_CPP) );
141 INTREPID_TEST_COMMAND( atools.contractFieldFieldScalar<
double>(a_10_2_2, a_10_2_2, a_10_2_3, COMP_CPP) );
142 INTREPID_TEST_COMMAND( atools.contractFieldFieldScalar<
double>(a_9_2_2, a_10_3_2, a_10_2_2, COMP_CPP) );
143 INTREPID_TEST_COMMAND( atools.contractFieldFieldScalar<
double>(a_10_3_2, a_10_2_2, a_10_3_2, COMP_CPP) );
144 INTREPID_TEST_COMMAND( atools.contractFieldFieldScalar<
double>(a_10_2_2, a_10_2_2, a_10_3_2, COMP_CPP) );
145 INTREPID_TEST_COMMAND( atools.contractFieldFieldScalar<
double>(a_10_2_3, a_10_2_2, a_10_3_2, COMP_ENGINE_MAX) );
146 INTREPID_TEST_COMMAND( atools.contractFieldFieldScalar<
double>(a_10_2_3, a_10_2_2, a_10_3_2, COMP_CPP) );
149 FieldContainer<double> a_10_2_2_2(10, 2, 2, 2);
150 FieldContainer<double> a_9_2_2_2(9, 2, 2, 2);
151 FieldContainer<double> a_10_3_2_2(10, 3, 2, 2);
152 FieldContainer<double> a_10_2_3_2(10, 2, 3, 2);
153 FieldContainer<double> a_10_2_2_3(10, 2, 2, 3);
156 *outStream <<
"-> contractFieldFieldVector:\n";
157 INTREPID_TEST_COMMAND( atools.contractFieldFieldVector<
double>(a_2_2, a_2_2, a_2_2, COMP_CPP) );
158 INTREPID_TEST_COMMAND( atools.contractFieldFieldVector<
double>(a_2_2, a_10_2_2_2, a_2_2, COMP_CPP) );
159 INTREPID_TEST_COMMAND( atools.contractFieldFieldVector<
double>(a_2_2, a_10_2_2_2, a_10_2_2_2, COMP_CPP) );
160 INTREPID_TEST_COMMAND( atools.contractFieldFieldVector<
double>(a_10_2_2, a_9_2_2_2, a_10_2_2_2, COMP_CPP) );
161 INTREPID_TEST_COMMAND( atools.contractFieldFieldVector<
double>(a_10_2_2, a_10_2_2_2, a_10_2_3_2, COMP_CPP) );
162 INTREPID_TEST_COMMAND( atools.contractFieldFieldVector<
double>(a_10_2_2, a_10_2_2_2, a_10_2_2_3, COMP_CPP) );
163 INTREPID_TEST_COMMAND( atools.contractFieldFieldVector<
double>(a_9_2_2, a_10_2_2_2, a_10_2_2_2, COMP_CPP) );
164 INTREPID_TEST_COMMAND( atools.contractFieldFieldVector<
double>(a_10_3_2, a_10_2_2_2, a_10_3_2_2, COMP_CPP) );
165 INTREPID_TEST_COMMAND( atools.contractFieldFieldVector<
double>(a_10_2_2, a_10_2_2_2, a_10_3_2_2, COMP_CPP) );
166 INTREPID_TEST_COMMAND( atools.contractFieldFieldVector<
double>(a_10_2_3, a_10_2_2_2, a_10_3_2_2, COMP_ENGINE_MAX) );
167 INTREPID_TEST_COMMAND( atools.contractFieldFieldVector<
double>(a_10_2_3, a_10_2_2_2, a_10_3_2_2, COMP_CPP) );
170 FieldContainer<double> a_10_2_2_2_2(10, 2, 2, 2, 2);
171 FieldContainer<double> a_9_2_2_2_2(9, 2, 2, 2, 2);
172 FieldContainer<double> a_10_3_2_2_2(10, 3, 2, 2, 2);
173 FieldContainer<double> a_10_2_3_2_2(10, 2, 3, 2, 2);
174 FieldContainer<double> a_10_2_2_3_2(10, 2, 2, 3, 2);
175 FieldContainer<double> a_10_2_2_2_3(10, 2, 2, 2, 3);
177 *outStream <<
"-> contractFieldFieldTensor:\n";
178 INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<
double>(a_2_2, a_2_2, a_2_2, COMP_CPP) );
179 INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<
double>(a_2_2, a_10_2_2_2_2, a_2_2, COMP_CPP) );
180 INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<
double>(a_2_2, a_10_2_2_2_2, a_10_2_2_2_2, COMP_CPP) );
181 INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<
double>(a_10_2_2, a_9_2_2_2_2, a_10_2_2_2_2, COMP_CPP) );
182 INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<
double>(a_10_2_2, a_10_2_2_2_2, a_10_2_3_2_2, COMP_CPP) );
183 INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<
double>(a_10_2_2, a_10_2_2_2_2, a_10_2_2_3_2, COMP_CPP) );
184 INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<
double>(a_10_2_2, a_10_2_2_2_2, a_10_2_2_2_3, COMP_CPP) );
185 INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<
double>(a_9_2_2, a_10_2_2_2_2, a_10_2_2_2_2, COMP_CPP) );
186 INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<
double>(a_10_3_2, a_10_2_2_2_2, a_10_2_2_2_2, COMP_CPP) );
187 INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<
double>(a_10_2_2, a_10_2_2_2_2, a_10_3_2_2_2, COMP_CPP) );
188 INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<
double>(a_10_2_3, a_10_2_2_2_2, a_10_3_2_2_2, COMP_ENGINE_MAX) );
189 INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<
double>(a_10_2_3, a_10_2_2_2_2, a_10_3_2_2_2, COMP_CPP) );
192 FieldContainer<double> a_9_2(9, 2);
193 FieldContainer<double> a_10_1(10, 1);
195 *outStream <<
"-> contractDataFieldScalar:\n";
196 INTREPID_TEST_COMMAND( atools.contractDataFieldScalar<
double>(a_2_2, a_2_2, a_2_2, COMP_CPP) );
197 INTREPID_TEST_COMMAND( atools.contractDataFieldScalar<
double>(a_2_2, a_10_2_2, a_10_2_2, COMP_CPP) );
198 INTREPID_TEST_COMMAND( atools.contractDataFieldScalar<
double>(a_10_2_2, a_2_2, a_10_2_2, COMP_CPP) );
199 INTREPID_TEST_COMMAND( atools.contractDataFieldScalar<
double>(a_2_2, a_10_2, a_9_2_2, COMP_CPP) );
200 INTREPID_TEST_COMMAND( atools.contractDataFieldScalar<
double>(a_2_2, a_10_3, a_10_2_2, COMP_CPP) );
201 INTREPID_TEST_COMMAND( atools.contractDataFieldScalar<
double>(a_9_2, a_10_2, a_10_2_2, COMP_CPP) );
202 INTREPID_TEST_COMMAND( atools.contractDataFieldScalar<
double>(a_10_2, a_10_2, a_10_3_2, COMP_CPP) );
203 INTREPID_TEST_COMMAND( atools.contractDataFieldScalar<
double>(a_10_2, a_10_2, a_10_2_2, COMP_ENGINE_MAX) );
204 INTREPID_TEST_COMMAND( atools.contractDataFieldScalar<
double>(a_10_2, a_10_2, a_10_2_2, COMP_CPP) );
205 INTREPID_TEST_COMMAND( atools.contractDataFieldScalar<
double>(a_10_2, a_10_1, a_10_2_2, COMP_CPP) );
208 FieldContainer<double> a_10_1_2(10, 1, 2);
209 FieldContainer<double> a_10_1_3(10, 1, 3);
211 *outStream <<
"-> contractDataFieldVector:\n";
212 INTREPID_TEST_COMMAND( atools.contractDataFieldVector<
double>(a_2_2, a_2_2, a_2_2, COMP_CPP) );
213 INTREPID_TEST_COMMAND( atools.contractDataFieldVector<
double>(a_2_2, a_2_2, a_10_2_2_2, COMP_CPP) );
214 INTREPID_TEST_COMMAND( atools.contractDataFieldVector<
double>(a_10_2_2, a_10_2_2, a_10_2_2_2, COMP_CPP) );
215 INTREPID_TEST_COMMAND( atools.contractDataFieldVector<
double>(a_2_2, a_10_2_2, a_9_2_2_2, COMP_CPP) );
216 INTREPID_TEST_COMMAND( atools.contractDataFieldVector<
double>(a_2_2, a_10_2_2, a_10_2_3_2, COMP_CPP) );
217 INTREPID_TEST_COMMAND( atools.contractDataFieldVector<
double>(a_2_2, a_10_2_2, a_10_2_2_3, COMP_CPP) );
218 INTREPID_TEST_COMMAND( atools.contractDataFieldVector<
double>(a_9_2, a_10_2_2, a_10_2_2_2, COMP_CPP) );
219 INTREPID_TEST_COMMAND( atools.contractDataFieldVector<
double>(a_10_2, a_10_2_2, a_10_3_2_2, COMP_CPP) );
220 INTREPID_TEST_COMMAND( atools.contractDataFieldVector<
double>(a_10_2, a_10_2_2, a_10_2_2_2, COMP_ENGINE_MAX) );
221 INTREPID_TEST_COMMAND( atools.contractDataFieldVector<
double>(a_10_2, a_10_2_2, a_10_2_2_2, COMP_CPP) );
222 INTREPID_TEST_COMMAND( atools.contractDataFieldVector<
double>(a_10_2, a_10_1_2, a_10_2_2_2, COMP_CPP) );
225 FieldContainer<double> a_10_1_2_2(10, 1, 2, 2);
227 *outStream <<
"-> contractDataFieldTensor:\n";
228 INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<
double>(a_2_2, a_2_2, a_2_2, COMP_CPP) );
229 INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<
double>(a_2_2, a_2_2, a_10_2_2_2_2, COMP_CPP) );
230 INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<
double>(a_10_2_2, a_10_2_2_2, a_10_2_2_2_2, COMP_CPP) );
231 INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<
double>(a_2_2, a_10_2_2_2, a_9_2_2_2_2, COMP_CPP) );
232 INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<
double>(a_2_2, a_10_2_2_2, a_10_2_3_2_2, COMP_CPP) );
233 INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<
double>(a_2_2, a_10_2_2_2, a_10_2_2_3_2, COMP_CPP) );
234 INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<
double>(a_2_2, a_10_2_2_2, a_10_2_2_2_3, COMP_CPP) );
235 INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<
double>(a_9_2, a_10_2_2_2, a_10_2_2_2_2, COMP_CPP) );
236 INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<
double>(a_10_2, a_10_2_2_2, a_10_3_2_2_2, COMP_CPP) );
237 INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<
double>(a_10_2, a_10_2_2_2, a_10_2_2_2_2, COMP_ENGINE_MAX) );
238 INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<
double>(a_10_2, a_10_2_2_2, a_10_2_2_2_2, COMP_CPP) );
239 INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<
double>(a_10_2, a_10_1_2_2, a_10_2_2_2_2, COMP_CPP) );
242 FieldContainer<double> a_2(2);
243 FieldContainer<double> a_10(10);
245 *outStream <<
"-> contractDataDataScalar:\n";
246 INTREPID_TEST_COMMAND( atools.contractDataDataScalar<
double>(a_2_2, a_10_2_2, a_10_2_2, COMP_CPP) );
247 INTREPID_TEST_COMMAND( atools.contractDataDataScalar<
double>(a_2_2, a_10_2, a_10_2_2, COMP_CPP) );
248 INTREPID_TEST_COMMAND( atools.contractDataDataScalar<
double>(a_2_2, a_10_2, a_10_2, COMP_CPP) );
249 INTREPID_TEST_COMMAND( atools.contractDataDataScalar<
double>(a_2, a_9_2, a_10_2, COMP_CPP) );
250 INTREPID_TEST_COMMAND( atools.contractDataDataScalar<
double>(a_2, a_10_2, a_10_3, COMP_CPP) );
251 INTREPID_TEST_COMMAND( atools.contractDataDataScalar<
double>(a_2, a_10_2, a_10_2, COMP_CPP) );
252 INTREPID_TEST_COMMAND( atools.contractDataDataScalar<
double>(a_10, a_10_2, a_10_2, COMP_ENGINE_MAX) );
253 INTREPID_TEST_COMMAND( atools.contractDataDataScalar<
double>(a_10, a_10_2, a_10_2, COMP_CPP) );
256 *outStream <<
"-> contractDataDataVector:\n";
257 INTREPID_TEST_COMMAND( atools.contractDataDataVector<
double>(a_2_2, a_2_2, a_2_2, COMP_CPP) );
258 INTREPID_TEST_COMMAND( atools.contractDataDataVector<
double>(a_2_2, a_10_2_2, a_2_2, COMP_CPP) );
259 INTREPID_TEST_COMMAND( atools.contractDataDataVector<
double>(a_2_2, a_10_2_2, a_10_2_2, COMP_CPP) );
260 INTREPID_TEST_COMMAND( atools.contractDataDataVector<
double>(a_10, a_9_2_2, a_10_2_2, COMP_CPP) );
261 INTREPID_TEST_COMMAND( atools.contractDataDataVector<
double>(a_10, a_10_3_2, a_10_2_2, COMP_CPP) );
262 INTREPID_TEST_COMMAND( atools.contractDataDataVector<
double>(a_10, a_10_2_3, a_10_2_2, COMP_CPP) );
263 INTREPID_TEST_COMMAND( atools.contractDataDataVector<
double>(a_2, a_10_2_2, a_10_2_2, COMP_CPP) );
264 INTREPID_TEST_COMMAND( atools.contractDataDataVector<
double>(a_10, a_10_2_2, a_10_2_2, COMP_ENGINE_MAX) );
265 INTREPID_TEST_COMMAND( atools.contractDataDataVector<
double>(a_10, a_10_2_2, a_10_2_2, COMP_CPP) );
268 *outStream <<
"-> contractDataDataTensor:\n";
269 INTREPID_TEST_COMMAND( atools.contractDataDataTensor<
double>(a_2_2, a_2_2, a_2_2, COMP_CPP) );
270 INTREPID_TEST_COMMAND( atools.contractDataDataTensor<
double>(a_2_2, a_10_2_2_2, a_2_2, COMP_CPP) );
271 INTREPID_TEST_COMMAND( atools.contractDataDataTensor<
double>(a_2_2, a_10_2_2_2, a_10_2_2_2, COMP_CPP) );
272 INTREPID_TEST_COMMAND( atools.contractDataDataTensor<
double>(a_10, a_9_2_2_2, a_10_2_2_2, COMP_CPP) );
273 INTREPID_TEST_COMMAND( atools.contractDataDataTensor<
double>(a_10, a_10_2_2_2, a_10_3_2_2, COMP_CPP) );
274 INTREPID_TEST_COMMAND( atools.contractDataDataTensor<
double>(a_10, a_10_2_2_2, a_10_2_3_2, COMP_CPP) );
275 INTREPID_TEST_COMMAND( atools.contractDataDataTensor<
double>(a_10, a_10_2_2_2, a_10_2_2_3, COMP_CPP) );
276 INTREPID_TEST_COMMAND( atools.contractDataDataTensor<
double>(a_2, a_10_2_2_2, a_10_2_2_2, COMP_CPP) );
277 INTREPID_TEST_COMMAND( atools.contractDataDataTensor<
double>(a_10, a_10_2_2_2, a_10_2_2_2, COMP_ENGINE_MAX) );
278 INTREPID_TEST_COMMAND( atools.contractDataDataTensor<
double>(a_10, a_10_2_2_2, a_10_2_2_2, COMP_CPP) );
284 catch (
const std::logic_error & err) {
285 *outStream <<
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
286 *outStream << err.what() <<
'\n';
287 *outStream <<
"-------------------------------------------------------------------------------" <<
"\n\n";
291#ifdef HAVE_INTREPID_DEBUG
292 if (Teuchos::TestForException_getThrowNumber() != endThrowNumber)
298 <<
"===============================================================================\n"\
299 <<
"| TEST 2: correctness of math operations |\n"\
300 <<
"===============================================================================\n";
302 outStream->precision(20);
306 *outStream <<
"\n************ Checking contractFieldFieldScalar ************\n";
308 int c=5, p=9, l=3, r=7;
310 FieldContainer<double> in_c_l_p(c, l, p);
311 FieldContainer<double> in_c_r_p(c, r, p);
312 FieldContainer<double> out1_c_l_r(c, l, r);
313 FieldContainer<double> out2_c_l_r(c, l, r);
314 double zero = INTREPID_TOL*10000.0;
317 for (
int i=0; i<in_c_l_p.size(); i++) {
318 in_c_l_p[i] = Teuchos::ScalarTraits<double>::random();
320 for (
int i=0; i<in_c_r_p.size(); i++) {
321 in_c_r_p[i] = Teuchos::ScalarTraits<double>::random();
324 art::contractFieldFieldScalar<double>(out1_c_l_r, in_c_l_p, in_c_r_p, COMP_CPP);
325 art::contractFieldFieldScalar<double>(out2_c_l_r, in_c_l_p, in_c_r_p, COMP_BLAS);
326 rst::subtract(&out1_c_l_r[0], &out2_c_l_r[0], out2_c_l_r.size());
327 if (rst::vectorNorm(&out1_c_l_r[0], out1_c_l_r.size(), NORM_ONE) > zero) {
328 *outStream <<
"\n\nINCORRECT contractFieldFieldScalar (1): check COMP_CPP vs. COMP_BLAS; "
329 <<
" diff-1norm = " << rst::vectorNorm(&out1_c_l_r[0], out1_c_l_r.size(), NORM_ONE) <<
"\n\n";
333 out1_c_l_r.initialize(2.0); out2_c_l_r.initialize(2.0);
334 art::contractFieldFieldScalar<double>(out1_c_l_r, in_c_l_p, in_c_r_p, COMP_CPP,
true);
335 art::contractFieldFieldScalar<double>(out2_c_l_r, in_c_l_p, in_c_r_p, COMP_BLAS,
true);
336 rst::subtract(&out1_c_l_r[0], &out2_c_l_r[0], out2_c_l_r.size());
337 if (rst::vectorNorm(&out1_c_l_r[0], out1_c_l_r.size(), NORM_ONE) > zero) {
338 *outStream <<
"\n\nINCORRECT contractFieldFieldScalar (1): check COMP_CPP vs. COMP_BLAS; "
339 <<
" diff-1norm = " << rst::vectorNorm(&out1_c_l_r[0], out1_c_l_r.size(), NORM_ONE) <<
"\n\n";
345 *outStream <<
"\n************ Checking contractFieldFieldVector ************\n";
347 int c=5, p=9, l=3, r=7, d=13;
349 FieldContainer<double> in_c_l_p_d(c, l, p, d);
350 FieldContainer<double> in_c_r_p_d(c, r, p, d);
351 FieldContainer<double> out1_c_l_r(c, l, r);
352 FieldContainer<double> out2_c_l_r(c, l, r);
353 double zero = INTREPID_TOL*10000.0;
356 for (
int i=0; i<in_c_l_p_d.size(); i++) {
357 in_c_l_p_d[i] = Teuchos::ScalarTraits<double>::random();
359 for (
int i=0; i<in_c_r_p_d.size(); i++) {
360 in_c_r_p_d[i] = Teuchos::ScalarTraits<double>::random();
363 art::contractFieldFieldVector<double>(out1_c_l_r, in_c_l_p_d, in_c_r_p_d, COMP_CPP);
364 art::contractFieldFieldVector<double>(out2_c_l_r, in_c_l_p_d, in_c_r_p_d, COMP_BLAS);
366 rst::subtract(&out1_c_l_r[0], &out2_c_l_r[0], out2_c_l_r.size());
367 if (rst::vectorNorm(&out1_c_l_r[0], out1_c_l_r.size(), NORM_ONE) > zero) {
368 *outStream <<
"\n\nINCORRECT contractFieldFieldVector (1): check COMP_CPP vs. COMP_BLAS; "
369 <<
" diff-1norm = " << rst::vectorNorm(&out1_c_l_r[0], out1_c_l_r.size(), NORM_ONE) <<
"\n\n";
374 out1_c_l_r.initialize(2.0); out2_c_l_r.initialize(2.0);
376 art::contractFieldFieldVector<double>(out1_c_l_r, in_c_l_p_d, in_c_r_p_d, COMP_CPP,
true);
377 art::contractFieldFieldVector<double>(out2_c_l_r, in_c_l_p_d, in_c_r_p_d, COMP_BLAS,
true);
379 rst::subtract(&out1_c_l_r[0], &out2_c_l_r[0], out2_c_l_r.size());
380 if (rst::vectorNorm(&out1_c_l_r[0], out1_c_l_r.size(), NORM_ONE) > zero) {
381 *outStream <<
"\n\nINCORRECT contractFieldFieldVector (1): check COMP_CPP vs. COMP_BLAS; "
382 <<
" diff-1norm = " << rst::vectorNorm(&out1_c_l_r[0], out1_c_l_r.size(), NORM_ONE) <<
"\n\n";
388 *outStream <<
"\n************ Checking contractFieldFieldTensor ************\n";
390 int c=5, p=9, l=3, r=7, d1=13, d2=5;
392 FieldContainer<double> in_c_l_p_d_d(c, l, p, d1, d2);
393 FieldContainer<double> in_c_r_p_d_d(c, r, p, d1, d2);
394 FieldContainer<double> out1_c_l_r(c, l, r);
395 FieldContainer<double> out2_c_l_r(c, l, r);
396 double zero = INTREPID_TOL*10000.0;
399 for (
int i=0; i<in_c_l_p_d_d.size(); i++) {
400 in_c_l_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
402 for (
int i=0; i<in_c_r_p_d_d.size(); i++) {
403 in_c_r_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
406 art::contractFieldFieldTensor<double>(out1_c_l_r, in_c_l_p_d_d, in_c_r_p_d_d, COMP_CPP);
407 art::contractFieldFieldTensor<double>(out2_c_l_r, in_c_l_p_d_d, in_c_r_p_d_d, COMP_BLAS);
409 rst::subtract(&out1_c_l_r[0], &out2_c_l_r[0], out2_c_l_r.size());
410 if (rst::vectorNorm(&out1_c_l_r[0], out1_c_l_r.size(), NORM_ONE) > zero) {
411 *outStream <<
"\n\nINCORRECT contractFieldFieldTensor (1): check COMP_CPP vs. COMP_BLAS; "
412 <<
" diff-1norm = " << rst::vectorNorm(&out1_c_l_r[0], out1_c_l_r.size(), NORM_ONE) <<
"\n\n";
417 out1_c_l_r.initialize(2.0); out2_c_l_r.initialize(2.0);
419 art::contractFieldFieldTensor<double>(out1_c_l_r, in_c_l_p_d_d, in_c_r_p_d_d, COMP_CPP,
true);
420 art::contractFieldFieldTensor<double>(out2_c_l_r, in_c_l_p_d_d, in_c_r_p_d_d, COMP_BLAS,
true);
422 rst::subtract(&out1_c_l_r[0], &out2_c_l_r[0], out2_c_l_r.size());
423 if (rst::vectorNorm(&out1_c_l_r[0], out1_c_l_r.size(), NORM_ONE) > zero) {
424 *outStream <<
"\n\nINCORRECT contractFieldFieldTensor (1): check COMP_CPP vs. COMP_BLAS; "
425 <<
" diff-1norm = " << rst::vectorNorm(&out1_c_l_r[0], out1_c_l_r.size(), NORM_ONE) <<
"\n\n";
431 *outStream <<
"\n************ Checking contractDataFieldScalar ************\n";
435 FieldContainer<double> in_c_l_p(c, l, p);
436 FieldContainer<double> data_c_p(c, p);
437 FieldContainer<double> data_c_1(c, 1);
438 FieldContainer<double> out1_c_l(c, l);
439 FieldContainer<double> out2_c_l(c, l);
440 double zero = INTREPID_TOL*10000.0;
443 for (
int i=0; i<in_c_l_p.size(); i++) {
444 in_c_l_p[i] = Teuchos::ScalarTraits<double>::random();
446 for (
int i=0; i<data_c_p.size(); i++) {
447 data_c_p[i] = Teuchos::ScalarTraits<double>::random();
449 for (
int i=0; i<data_c_1.size(); i++) {
450 data_c_1[i] = Teuchos::ScalarTraits<double>::random();
454 art::contractDataFieldScalar<double>(out1_c_l, data_c_p, in_c_l_p, COMP_CPP);
455 art::contractDataFieldScalar<double>(out2_c_l, data_c_p, in_c_l_p, COMP_BLAS);
456 rst::subtract(&out1_c_l[0], &out2_c_l[0], out2_c_l.size());
457 if (rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) > zero) {
458 *outStream <<
"\n\nINCORRECT contractDataFieldScalar (1): check COMP_CPP vs. COMP_BLAS; "
459 <<
" diff-1norm = " << rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) <<
"\n\n";
463 art::contractDataFieldScalar<double>(out1_c_l, data_c_1, in_c_l_p, COMP_CPP);
464 art::contractDataFieldScalar<double>(out2_c_l, data_c_1, in_c_l_p, COMP_BLAS);
465 rst::subtract(&out1_c_l[0], &out2_c_l[0], out2_c_l.size());
466 if (rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) > zero) {
467 *outStream <<
"\n\nINCORRECT contractDataFieldScalar (2): check COMP_CPP vs. COMP_BLAS; "
468 <<
" diff-1norm = " << rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) <<
"\n\n";
472 out1_c_l.initialize(2.0); out2_c_l.initialize(2.0);
473 art::contractDataFieldScalar<double>(out1_c_l, data_c_p, in_c_l_p, COMP_CPP,
true);
474 art::contractDataFieldScalar<double>(out2_c_l, data_c_p, in_c_l_p, COMP_BLAS,
true);
475 rst::subtract(&out1_c_l[0], &out2_c_l[0], out2_c_l.size());
476 if (rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) > zero) {
477 *outStream <<
"\n\nINCORRECT contractDataFieldScalar (1): check COMP_CPP vs. COMP_BLAS; "
478 <<
" diff-1norm = " << rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) <<
"\n\n";
484 *outStream <<
"\n************ Checking contractDataFieldVector ************\n";
486 int c=5, p=9, l=7, d=3;
488 FieldContainer<double> in_c_l_p_d(c, l, p, d);
489 FieldContainer<double> data_c_p_d(c, p, d);
490 FieldContainer<double> data_c_1_d(c, 1, d);
491 FieldContainer<double> out1_c_l(c, l);
492 FieldContainer<double> out2_c_l(c, l);
493 double zero = INTREPID_TOL*10000.0;
496 for (
int i=0; i<in_c_l_p_d.size(); i++) {
497 in_c_l_p_d[i] = Teuchos::ScalarTraits<double>::random();
499 for (
int i=0; i<data_c_p_d.size(); i++) {
500 data_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
502 for (
int i=0; i<data_c_1_d.size(); i++) {
503 data_c_1_d[i] = Teuchos::ScalarTraits<double>::random();
507 art::contractDataFieldVector<double>(out1_c_l, data_c_p_d, in_c_l_p_d, COMP_CPP);
508 art::contractDataFieldVector<double>(out2_c_l, data_c_p_d, in_c_l_p_d, COMP_BLAS);
509 rst::subtract(&out1_c_l[0], &out2_c_l[0], out2_c_l.size());
510 if (rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) > zero) {
511 *outStream <<
"\n\nINCORRECT contractDataFieldVector (1): check COMP_CPP vs. COMP_BLAS; "
512 <<
" diff-1norm = " << rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) <<
"\n\n";
516 art::contractDataFieldVector<double>(out1_c_l, data_c_1_d, in_c_l_p_d, COMP_CPP);
517 art::contractDataFieldVector<double>(out2_c_l, data_c_1_d, in_c_l_p_d, COMP_BLAS);
518 rst::subtract(&out1_c_l[0], &out2_c_l[0], out2_c_l.size());
519 if (rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) > zero) {
520 *outStream <<
"\n\nINCORRECT contractDataFieldVector (2): check COMP_CPP vs. COMP_BLAS; "
521 <<
" diff-1norm = " << rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) <<
"\n\n";
525 out1_c_l.initialize(2.0); out2_c_l.initialize(2.0);
526 art::contractDataFieldVector<double>(out1_c_l, data_c_p_d, in_c_l_p_d, COMP_CPP,
true);
527 art::contractDataFieldVector<double>(out2_c_l, data_c_p_d, in_c_l_p_d, COMP_BLAS,
true);
528 rst::subtract(&out1_c_l[0], &out2_c_l[0], out2_c_l.size());
529 if (rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) > zero) {
530 *outStream <<
"\n\nINCORRECT contractDataFieldVector (1): check COMP_CPP vs. COMP_BLAS; "
531 <<
" diff-1norm = " << rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) <<
"\n\n";
537 *outStream <<
"\n************ Checking contractDataFieldTensor ************\n";
539 int c=5, p=9, l=7, d1=3, d2=13;
541 FieldContainer<double> in_c_l_p_d_d(c, l, p, d1, d2);
542 FieldContainer<double> data_c_p_d_d(c, p, d1, d2);
543 FieldContainer<double> data_c_1_d_d(c, 1, d1, d2);
544 FieldContainer<double> out1_c_l(c, l);
545 FieldContainer<double> out2_c_l(c, l);
546 double zero = INTREPID_TOL*10000.0;
549 for (
int i=0; i<in_c_l_p_d_d.size(); i++) {
550 in_c_l_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
552 for (
int i=0; i<data_c_p_d_d.size(); i++) {
553 data_c_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
555 for (
int i=0; i<data_c_1_d_d.size(); i++) {
556 data_c_1_d_d[i] = Teuchos::ScalarTraits<double>::random();
560 art::contractDataFieldTensor<double>(out1_c_l, data_c_p_d_d, in_c_l_p_d_d, COMP_CPP);
561 art::contractDataFieldTensor<double>(out2_c_l, data_c_p_d_d, in_c_l_p_d_d, COMP_BLAS);
562 rst::subtract(&out1_c_l[0], &out2_c_l[0], out2_c_l.size());
563 if (rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) > zero) {
564 *outStream <<
"\n\nINCORRECT contractDataFieldTensor (1): check COMP_CPP vs. COMP_BLAS; "
565 <<
" diff-1norm = " << rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) <<
"\n\n";
569 art::contractDataFieldTensor<double>(out1_c_l, data_c_1_d_d, in_c_l_p_d_d, COMP_CPP);
570 art::contractDataFieldTensor<double>(out2_c_l, data_c_1_d_d, in_c_l_p_d_d, COMP_BLAS);
571 rst::subtract(&out1_c_l[0], &out2_c_l[0], out2_c_l.size());
572 if (rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) > zero) {
573 *outStream <<
"\n\nINCORRECT contractDataFieldTensor (1): check COMP_CPP vs. COMP_BLAS; "
574 <<
" diff-1norm = " << rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) <<
"\n\n";
578 out1_c_l.initialize(2.0); out2_c_l.initialize(2.0);
579 art::contractDataFieldTensor<double>(out1_c_l, data_c_p_d_d, in_c_l_p_d_d, COMP_CPP,
true);
580 art::contractDataFieldTensor<double>(out2_c_l, data_c_p_d_d, in_c_l_p_d_d, COMP_BLAS,
true);
581 rst::subtract(&out1_c_l[0], &out2_c_l[0], out2_c_l.size());
582 if (rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) > zero) {
583 *outStream <<
"\n\nINCORRECT contractDataFieldTensor (1): check COMP_CPP vs. COMP_BLAS; "
584 <<
" diff-1norm = " << rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) <<
"\n\n";
590 *outStream <<
"\n************ Checking contractDataDataScalar ************\n";
594 FieldContainer<double> inl_c_p(c, p);
595 FieldContainer<double> inr_c_p(c, p);
596 FieldContainer<double> out1_c(c);
597 FieldContainer<double> out2_c(c);
598 double zero = INTREPID_TOL*10000.0;
601 for (
int i=0; i<inl_c_p.size(); i++) {
602 inl_c_p[i] = Teuchos::ScalarTraits<double>::random();
604 for (
int i=0; i<inr_c_p.size(); i++) {
605 inr_c_p[i] = Teuchos::ScalarTraits<double>::random();
608 art::contractDataDataScalar<double>(out1_c, inl_c_p, inr_c_p, COMP_CPP);
609 art::contractDataDataScalar<double>(out2_c, inl_c_p, inr_c_p, COMP_BLAS);
610 rst::subtract(&out1_c[0], &out2_c[0], out2_c.size());
611 if (rst::vectorNorm(&out1_c[0], out1_c.size(), NORM_ONE) > zero) {
612 *outStream <<
"\n\nINCORRECT contractDataDataScalar (1): check COMP_CPP vs. COMP_BLAS; "
613 <<
" diff-1norm = " << rst::vectorNorm(&out1_c[0], out1_c.size(), NORM_ONE) <<
"\n\n";
617 out1_c.initialize(2.0); out2_c.initialize(2.0);
618 art::contractDataDataScalar<double>(out1_c, inl_c_p, inr_c_p, COMP_CPP,
true);
619 art::contractDataDataScalar<double>(out2_c, inl_c_p, inr_c_p, COMP_BLAS,
true);
620 rst::subtract(&out1_c[0], &out2_c[0], out2_c.size());
621 if (rst::vectorNorm(&out1_c[0], out1_c.size(), NORM_ONE) > zero) {
622 *outStream <<
"\n\nINCORRECT contractDataDataScalar (1): check COMP_CPP vs. COMP_BLAS; "
623 <<
" diff-1norm = " << rst::vectorNorm(&out1_c[0], out1_c.size(), NORM_ONE) <<
"\n\n";
629 *outStream <<
"\n************ Checking contractDataDataVector ************\n";
633 FieldContainer<double> inl_c_p_d(c, p, d);
634 FieldContainer<double> inr_c_p_d(c, p, d);
635 FieldContainer<double> out1_c(c);
636 FieldContainer<double> out2_c(c);
637 double zero = INTREPID_TOL*10000.0;
640 for (
int i=0; i<inl_c_p_d.size(); i++) {
641 inl_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
643 for (
int i=0; i<inr_c_p_d.size(); i++) {
644 inr_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
647 art::contractDataDataVector<double>(out1_c, inl_c_p_d, inr_c_p_d, COMP_CPP);
648 art::contractDataDataVector<double>(out2_c, inl_c_p_d, inr_c_p_d, COMP_BLAS);
650 rst::subtract(&out1_c[0], &out2_c[0], out2_c.size());
651 if (rst::vectorNorm(&out1_c[0], out1_c.size(), NORM_ONE) > zero) {
652 *outStream <<
"\n\nINCORRECT contractDataDataVector (1): check COMP_CPP vs. COMP_BLAS; "
653 <<
" diff-1norm = " << rst::vectorNorm(&out1_c[0], out1_c.size(), NORM_ONE) <<
"\n\n";
658 out1_c.initialize(2.0); out2_c.initialize(2.0);
660 art::contractDataDataVector<double>(out1_c, inl_c_p_d, inr_c_p_d, COMP_CPP,
true);
661 art::contractDataDataVector<double>(out2_c, inl_c_p_d, inr_c_p_d, COMP_BLAS,
true);
663 rst::subtract(&out1_c[0], &out2_c[0], out2_c.size());
664 if (rst::vectorNorm(&out1_c[0], out1_c.size(), NORM_ONE) > zero) {
665 *outStream <<
"\n\nINCORRECT contractDataDataVector (1): check COMP_CPP vs. COMP_BLAS; "
666 <<
" diff-1norm = " << rst::vectorNorm(&out1_c[0], out1_c.size(), NORM_ONE) <<
"\n\n";
672 *outStream <<
"\n************ Checking contractDataDataTensor ************\n";
674 int c=5, p=9, d1=13, d2=5;
676 FieldContainer<double> inl_c_p_d_d(c, p, d1, d2);
677 FieldContainer<double> inr_c_p_d_d(c, p, d1, d2);
678 FieldContainer<double> out1_c(c);
679 FieldContainer<double> out2_c(c);
680 double zero = INTREPID_TOL*10000.0;
683 for (
int i=0; i<inl_c_p_d_d.size(); i++) {
684 inl_c_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
686 for (
int i=0; i<inr_c_p_d_d.size(); i++) {
687 inr_c_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
690 art::contractDataDataTensor<double>(out1_c, inl_c_p_d_d, inr_c_p_d_d, COMP_CPP);
691 art::contractDataDataTensor<double>(out2_c, inl_c_p_d_d, inr_c_p_d_d, COMP_BLAS);
693 rst::subtract(&out1_c[0], &out2_c[0], out2_c.size());
694 if (rst::vectorNorm(&out1_c[0], out1_c.size(), NORM_ONE) > zero) {
695 *outStream <<
"\n\nINCORRECT contractDataDataTensor (1): check COMP_CPP vs. COMP_BLAS; "
696 <<
" diff-1norm = " << rst::vectorNorm(&out1_c[0], out1_c.size(), NORM_ONE) <<
"\n\n";
701 out1_c.initialize(2.0); out2_c.initialize(2.0);
703 art::contractDataDataTensor<double>(out1_c, inl_c_p_d_d, inr_c_p_d_d, COMP_CPP,
true);
704 art::contractDataDataTensor<double>(out2_c, inl_c_p_d_d, inr_c_p_d_d, COMP_BLAS,
true);
706 rst::subtract(&out1_c[0], &out2_c[0], out2_c.size());
707 if (rst::vectorNorm(&out1_c[0], out1_c.size(), NORM_ONE) > zero) {
708 *outStream <<
"\n\nINCORRECT contractDataDataTensor (1): check COMP_CPP vs. COMP_BLAS; "
709 <<
" diff-1norm = " << rst::vectorNorm(&out1_c[0], out1_c.size(), NORM_ONE) <<
"\n\n";
717 catch (
const std::logic_error & err) {
718 *outStream <<
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
719 *outStream << err.what() <<
'\n';
720 *outStream <<
"-------------------------------------------------------------------------------" <<
"\n\n";
726 std::cout <<
"End Result: TEST FAILED\n";
728 std::cout <<
"End Result: TEST PASSED\n";
731 std::cout.copyfmt(oldFormatState);
Header file for utility class to provide multidimensional containers.