Intrepid
test_01.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
49#include "Intrepid_RealSpaceTools.hpp"
50#include "Intrepid_FieldContainer.hpp"
51#include "Teuchos_oblackholestream.hpp"
52#include "Teuchos_RCP.hpp"
53#include "Teuchos_ScalarTraits.hpp"
54#include "Teuchos_GlobalMPISession.hpp"
55
56
57using namespace std;
58using namespace Intrepid;
59
60int main(int argc, char *argv[]) {
61
62 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
63
64 // This little trick lets us print to std::cout only if
65 // a (dummy) command-line argument is provided.
66 int iprint = argc - 1;
67 Teuchos::RCP<std::ostream> outStream;
68 Teuchos::oblackholestream bhs; // outputs nothing
69 if (iprint > 0)
70 outStream = Teuchos::rcp(&std::cout, false);
71 else
72 outStream = Teuchos::rcp(&bhs, false);
73
74 // Save the format state of the original std::cout.
75 Teuchos::oblackholestream oldFormatState;
76 oldFormatState.copyfmt(std::cout);
77
78 *outStream \
79 << "===============================================================================\n" \
80 << "| |\n" \
81 << "| Unit Test (RealSpaceTools) |\n" \
82 << "| |\n" \
83 << "| 1) Vector operations in 1D, 2D, and 3D real space |\n" \
84 << "| 2) Matrix / matrix-vector operations in 1D, 2D, and 3D real space |\n" \
85 << "| |\n" \
86 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \
87 << "| Denis Ridzal (dridzal@sandia.gov). |\n" \
88 << "| |\n" \
89 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
90 << "| Trilinos website: http://trilinos.sandia.gov |\n" \
91 << "| |\n" \
92 << "===============================================================================\n";
93
94
95 typedef RealSpaceTools<double> rst;
96
97
98 int errorFlag = 0;
99#ifdef HAVE_INTREPID_DEBUG
100 int beginThrowNumber = Teuchos::TestForException_getThrowNumber();
101 int endThrowNumber = beginThrowNumber + 49;
102#ifndef HAVE_INTREPID_DEBUG_INF_CHECK
103 endThrowNumber = beginThrowNumber + 44;
104#endif
105
106#endif
107
108 *outStream \
109 << "\n"
110 << "===============================================================================\n"\
111 << "| TEST 1: vector exceptions |\n"\
112 << "===============================================================================\n";
113
114 try{
115
116 FieldContainer<double> a_2_2(2, 2);
117 FieldContainer<double> a_10_2(10, 2);
118 FieldContainer<double> a_10_3(10, 3);
119 FieldContainer<double> a_10_2_2(10, 2, 2);
120 FieldContainer<double> a_10_2_3(10, 2, 3);
121 FieldContainer<double> a_10_2_2_3(10, 2, 2, 3);
122
123#ifdef HAVE_INTREPID_DEBUG
124
125 *outStream << "-> vector norm with multidimensional arrays:\n";
126
127 try {
128 rst::vectorNorm(a_2_2, NORM_TWO);
129 }
130 catch (const std::logic_error & err) {
131 *outStream << "Expected Error ----------------------------------------------------------------\n";
132 *outStream << err.what() << '\n';
133 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
134 };
135 try {
136 rst::vectorNorm(a_10_2_2, a_10_2_2, NORM_TWO);
137 }
138 catch (const std::logic_error & err) {
139 *outStream << "Expected Error ----------------------------------------------------------------\n";
140 *outStream << err.what() << '\n';
141 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
142 };
143 try {
144 rst::vectorNorm(a_10_2_2, a_10_2_2_3, NORM_TWO);
145 }
146 catch (const std::logic_error & err) {
147 *outStream << "Expected Error ----------------------------------------------------------------\n";
148 *outStream << err.what() << '\n';
149 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
150 };
151 try {
152 rst::vectorNorm(a_10_3, a_10_2_2, NORM_TWO);
153 }
154 catch (const std::logic_error & err) {
155 *outStream << "Expected Error ----------------------------------------------------------------\n";
156 *outStream << err.what() << '\n';
157 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
158 };
159
160 *outStream << "-> add with multidimensional arrays:\n";
161
162 try {
163 rst::add(a_10_2_2, a_10_2, a_2_2);
164 }
165 catch (const std::logic_error & err) {
166 *outStream << "Expected Error ----------------------------------------------------------------\n";
167 *outStream << err.what() << '\n';
168 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
169 };
170 try {
171 rst::add(a_10_2_3, a_10_2_2, a_10_2_2);
172 }
173 catch (const std::logic_error & err) {
174 *outStream << "Expected Error ----------------------------------------------------------------\n";
175 *outStream << err.what() << '\n';
176 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
177 };
178 try {
179 rst::add(a_10_2_2, a_10_2_2_3);
180 }
181 catch (const std::logic_error & err) {
182 *outStream << "Expected Error ----------------------------------------------------------------\n";
183 *outStream << err.what() << '\n';
184 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
185 };
186 try {
187 rst::add(a_10_2_3, a_10_2_2);
188 }
189 catch (const std::logic_error & err) {
190 *outStream << "Expected Error ----------------------------------------------------------------\n";
191 *outStream << err.what() << '\n';
192 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
193 };
194
195 *outStream << "-> subtract with multidimensional arrays:\n";
196
197 try {
198 rst::subtract(a_10_2_2, a_10_2, a_2_2);
199 }
200 catch (const std::logic_error & err) {
201 *outStream << "Expected Error ----------------------------------------------------------------\n";
202 *outStream << err.what() << '\n';
203 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
204 };
205 try {
206 rst::subtract(a_10_2_3, a_10_2_2, a_10_2_2);
207 }
208 catch (const std::logic_error & err) {
209 *outStream << "Expected Error ----------------------------------------------------------------\n";
210 *outStream << err.what() << '\n';
211 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
212 };
213 try {
214 rst::subtract(a_10_2_2, a_10_2_2_3);
215 }
216 catch (const std::logic_error & err) {
217 *outStream << "Expected Error ----------------------------------------------------------------\n";
218 *outStream << err.what() << '\n';
219 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
220 };
221 try {
222 rst::subtract(a_10_2_3, a_10_2_2);
223 }
224 catch (const std::logic_error & err) {
225 *outStream << "Expected Error ----------------------------------------------------------------\n";
226 *outStream << err.what() << '\n';
227 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
228 };
229
230 *outStream << "-> dot product norm with multidimensional arrays:\n";
231
232 try {
233 rst::dot(a_10_2, a_10_2_2_3, a_10_2_2_3);
234 }
235 catch (const std::logic_error & err) {
236 *outStream << "Expected Error ----------------------------------------------------------------\n";
237 *outStream << err.what() << '\n';
238 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
239 };
240 try {
241 rst::dot(a_10_2, a_10_2_2, a_10_2_2_3);
242 }
243 catch (const std::logic_error & err) {
244 *outStream << "Expected Error ----------------------------------------------------------------\n";
245 *outStream << err.what() << '\n';
246 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
247 };
248 try {
249 rst::dot(a_10_2_2, a_10_2_2_3, a_10_2_2_3);
250 }
251 catch (const std::logic_error & err) {
252 *outStream << "Expected Error ----------------------------------------------------------------\n";
253 *outStream << err.what() << '\n';
254 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
255 };
256 try {
257 rst::dot(a_10_2, a_10_2_2, a_10_2_3);
258 }
259 catch (const std::logic_error & err) {
260 *outStream << "Expected Error ----------------------------------------------------------------\n";
261 *outStream << err.what() << '\n';
262 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
263 };
264 try {
265 rst::dot(a_10_3, a_10_2_3, a_10_2_3);
266 }
267 catch (const std::logic_error & err) {
268 *outStream << "Expected Error ----------------------------------------------------------------\n";
269 *outStream << err.what() << '\n';
270 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
271 };
272
273 *outStream << "-> absolute value with multidimensional arrays:\n";
274
275 try {
276 rst::absval(a_10_3, a_10_2_3);
277 }
278 catch (const std::logic_error & err) {
279 *outStream << "Expected Error ----------------------------------------------------------------\n";
280 *outStream << err.what() << '\n';
281 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
282 };
283 try {
284 rst::absval(a_10_2_2, a_10_2_3);
285 }
286 catch (const std::logic_error & err) {
287 *outStream << "Expected Error ----------------------------------------------------------------\n";
288 *outStream << err.what() << '\n';
289 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
290 };
291#endif
292
293 }
294 catch (const std::logic_error & err) {
295 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
296 *outStream << err.what() << '\n';
297 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
298 errorFlag = -1000;
299 };
300
301
302
303 *outStream \
304 << "\n"
305 << "===============================================================================\n"\
306 << "| TEST 2: matrix / matrix-vector exceptions |\n"\
307 << "===============================================================================\n";
308
309 try{
310
311 FieldContainer<double> a_10_1_2_3_4(10, 1, 2, 3, 4);
312 FieldContainer<double> b_10_1_2_3_4(10, 1, 2, 3, 4);
313 FieldContainer<double> a_10(10);
314 FieldContainer<double> a_9(9);
315 FieldContainer<double> b_10(10);
316 FieldContainer<double> a_10_15_4_4(10, 15, 4, 4);
317 FieldContainer<double> b_10_15_4_4(10, 15, 4, 4);
318 FieldContainer<double> a_10_2_2(10, 2, 2);
319 FieldContainer<double> a_10_2_3(10, 2, 3);
320 FieldContainer<double> b_10_2_3(10, 2, 3);
321
322 FieldContainer<double> a_1_1(1, 1);
323 FieldContainer<double> b_1_1(1, 1);
324 FieldContainer<double> a_2_2(2, 2);
325 FieldContainer<double> b_2_2(2, 2);
326 FieldContainer<double> a_3_3(3, 3);
327 FieldContainer<double> b_3_3(3, 3);
328 FieldContainer<double> a_2_3(2, 3);
329 FieldContainer<double> a_4_4(4, 4);
330
331 FieldContainer<double> a_10_15_1_1(10, 15, 1, 1);
332 FieldContainer<double> b_10_15_1_1(10, 15, 1, 1);
333 FieldContainer<double> a_10_15_2_2(10, 15, 2, 2);
334 FieldContainer<double> b_10_15_2_2(10, 15, 2, 2);
335 FieldContainer<double> a_10_15_3_3(10, 15, 3, 3);
336 FieldContainer<double> a_10_15_3_2(10, 15, 3, 2);
337 FieldContainer<double> b_10_15_3_3(10, 15, 3, 3);
338 FieldContainer<double> b_10_14(10, 14);
339 FieldContainer<double> b_10_15(10, 15);
340 FieldContainer<double> b_10_14_3(10, 14, 3);
341 FieldContainer<double> b_10_15_3(10, 15, 3);
342
343
344#ifdef HAVE_INTREPID_DEBUG
345
346 *outStream << "-> inverse with multidimensional arrays:\n";
347
348 try {
349 rst::inverse(a_2_2, a_10_2_2);
350 }
351 catch (const std::logic_error & err) {
352 *outStream << "Expected Error ----------------------------------------------------------------\n";
353 *outStream << err.what() << '\n';
354 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
355 };
356 try {
357 rst::inverse(b_10_1_2_3_4, a_10_1_2_3_4);
358 }
359 catch (const std::logic_error & err) {
360 *outStream << "Expected Error ----------------------------------------------------------------\n";
361 *outStream << err.what() << '\n';
362 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
363 };
364 try {
365 rst::inverse(b_10, a_10);
366 }
367 catch (const std::logic_error & err) {
368 *outStream << "Expected Error ----------------------------------------------------------------\n";
369 *outStream << err.what() << '\n';
370 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
371 };
372 try {
373 rst::inverse(a_10_2_2, a_10_2_3);
374 }
375 catch (const std::logic_error & err) {
376 *outStream << "Expected Error ----------------------------------------------------------------\n";
377 *outStream << err.what() << '\n';
378 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
379 };
380 try {
381 rst::inverse(b_10_2_3, a_10_2_3);
382 }
383 catch (const std::logic_error & err) {
384 *outStream << "Expected Error ----------------------------------------------------------------\n";
385 *outStream << err.what() << '\n';
386 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
387 };
388 try {
389 rst::inverse(b_10_15_4_4, a_10_15_4_4);
390 }
391 catch (const std::logic_error & err) {
392 *outStream << "Expected Error ----------------------------------------------------------------\n";
393 *outStream << err.what() << '\n';
394 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
395 };
396 try {
397 rst::inverse(b_1_1, a_1_1);
398 }
399 catch (const std::logic_error & err) {
400 *outStream << "Expected Error ----------------------------------------------------------------\n";
401 *outStream << err.what() << '\n';
402 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
403 };
404 try {
405 rst::inverse(b_2_2, a_2_2);
406 }
407 catch (const std::logic_error & err) {
408 *outStream << "Expected Error ----------------------------------------------------------------\n";
409 *outStream << err.what() << '\n';
410 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
411 };
412 try {
413 rst::inverse(b_3_3, a_3_3);
414 }
415 catch (const std::logic_error & err) {
416 *outStream << "Expected Error ----------------------------------------------------------------\n";
417 *outStream << err.what() << '\n';
418 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
419 };
420 a_2_2[0] = 1.0;
421 a_3_3[0] = 1.0;
422 try {
423 rst::inverse(b_2_2, a_2_2);
424 }
425 catch (const std::logic_error & err) {
426 *outStream << "Expected Error ----------------------------------------------------------------\n";
427 *outStream << err.what() << '\n';
428 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
429 };
430 try {
431 rst::inverse(b_3_3, a_3_3);
432 }
433 catch (const std::logic_error & err) {
434 *outStream << "Expected Error ----------------------------------------------------------------\n";
435 *outStream << err.what() << '\n';
436 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
437 };
438
439 *outStream << "-> transpose with multidimensional arrays:\n";
440
441 try {
442 rst::transpose(a_2_2, a_10_2_2);
443 }
444 catch (const std::logic_error & err) {
445 *outStream << "Expected Error ----------------------------------------------------------------\n";
446 *outStream << err.what() << '\n';
447 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
448 };
449 try {
450 rst::transpose(b_10_1_2_3_4, a_10_1_2_3_4);
451 }
452 catch (const std::logic_error & err) {
453 *outStream << "Expected Error ----------------------------------------------------------------\n";
454 *outStream << err.what() << '\n';
455 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
456 };
457 try {
458 rst::transpose(b_10, a_10);
459 }
460 catch (const std::logic_error & err) {
461 *outStream << "Expected Error ----------------------------------------------------------------\n";
462 *outStream << err.what() << '\n';
463 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
464 };
465 try {
466 rst::transpose(a_10_2_2, a_10_2_3);
467 }
468 catch (const std::logic_error & err) {
469 *outStream << "Expected Error ----------------------------------------------------------------\n";
470 *outStream << err.what() << '\n';
471 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
472 };
473 try {
474 rst::transpose(b_10_2_3, a_10_2_3);
475 }
476 catch (const std::logic_error & err) {
477 *outStream << "Expected Error ----------------------------------------------------------------\n";
478 *outStream << err.what() << '\n';
479 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
480 };
481
482 *outStream << "-> determinant with multidimensional arrays:\n";
483
484 try {
485 rst::det(a_2_2, a_10_2_2);
486 }
487 catch (const std::logic_error & err) {
488 *outStream << "Expected Error ----------------------------------------------------------------\n";
489 *outStream << err.what() << '\n';
490 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
491 };
492 try {
493 rst::det(a_10_2_2, a_10_1_2_3_4);
494 }
495 catch (const std::logic_error & err) {
496 *outStream << "Expected Error ----------------------------------------------------------------\n";
497 *outStream << err.what() << '\n';
498 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
499 };
500 try {
501 rst::det(b_10_14, a_10_15_3_3);
502 }
503 catch (const std::logic_error & err) {
504 *outStream << "Expected Error ----------------------------------------------------------------\n";
505 *outStream << err.what() << '\n';
506 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
507 };
508 try {
509 rst::det(a_9, a_10_2_2);
510 }
511 catch (const std::logic_error & err) {
512 *outStream << "Expected Error ----------------------------------------------------------------\n";
513 *outStream << err.what() << '\n';
514 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
515 };
516 try {
517 rst::det(b_10, a_10_2_3);
518 }
519 catch (const std::logic_error & err) {
520 *outStream << "Expected Error ----------------------------------------------------------------\n";
521 *outStream << err.what() << '\n';
522 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
523 };
524 try {
525 rst::det(b_10_15, a_10_15_4_4);
526 }
527 catch (const std::logic_error & err) {
528 *outStream << "Expected Error ----------------------------------------------------------------\n";
529 *outStream << err.what() << '\n';
530 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
531 };
532 try {
533 rst::det(a_10_15_4_4);
534 }
535 catch (const std::logic_error & err) {
536 *outStream << "Expected Error ----------------------------------------------------------------\n";
537 *outStream << err.what() << '\n';
538 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
539 };
540 try {
541 rst::det(a_2_3);
542 }
543 catch (const std::logic_error & err) {
544 *outStream << "Expected Error ----------------------------------------------------------------\n";
545 *outStream << err.what() << '\n';
546 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
547 };
548 try {
549 rst::det(a_4_4);
550 }
551 catch (const std::logic_error & err) {
552 *outStream << "Expected Error ----------------------------------------------------------------\n";
553 *outStream << err.what() << '\n';
554 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
555 };
556
557 *outStream << "-> matrix-vector product with multidimensional arrays:\n";
558
559 try {
560 rst::matvec(a_10_2_2, a_10_2_3, b_10_2_3);
561 }
562 catch (const std::logic_error & err) {
563 *outStream << "Expected Error ----------------------------------------------------------------\n";
564 *outStream << err.what() << '\n';
565 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
566 };
567 try {
568 rst::matvec(a_2_2, a_2_2, a_10);
569 }
570 catch (const std::logic_error & err) {
571 *outStream << "Expected Error ----------------------------------------------------------------\n";
572 *outStream << err.what() << '\n';
573 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
574 };
575 try {
576 rst::matvec(a_9, a_10_2_2, a_2_2);
577 }
578 catch (const std::logic_error & err) {
579 *outStream << "Expected Error ----------------------------------------------------------------\n";
580 *outStream << err.what() << '\n';
581 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
582 };
583 try {
584 rst::matvec(b_10_15_3, a_10_15_3_3, b_10_14_3);
585 }
586 catch (const std::logic_error & err) {
587 *outStream << "Expected Error ----------------------------------------------------------------\n";
588 *outStream << err.what() << '\n';
589 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
590 };
591 try {
592 rst::matvec(b_10_14_3, a_10_15_3_3, b_10_15_3);
593 }
594 catch (const std::logic_error & err) {
595 *outStream << "Expected Error ----------------------------------------------------------------\n";
596 *outStream << err.what() << '\n';
597 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
598 };
599 try {
600 rst::matvec(b_10_15_3, a_10_15_3_2, b_10_15_3);
601 }
602 catch (const std::logic_error & err) {
603 *outStream << "Expected Error ----------------------------------------------------------------\n";
604 *outStream << err.what() << '\n';
605 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
606 };
607
608#endif
609
610 }
611 catch (const std::logic_error & err) {
612 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
613 *outStream << err.what() << '\n';
614 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
615 errorFlag = -1000;
616 };
617#ifdef HAVE_INTREPID_DEBUG
618 if (Teuchos::TestForException_getThrowNumber() != endThrowNumber)
619 errorFlag++;
620#endif
621
622
623 *outStream \
624 << "\n"
625 << "===============================================================================\n"\
626 << "| TEST 2: correctness of math operations |\n"\
627 << "===============================================================================\n";
628
629 outStream->precision(20);
630
631 try{
632
633 // two-dimensional base containers
634 for (int dim=3; dim>0; dim--) {
635 int i0=4, i1=5;
636 FieldContainer<double> ma_x_x_d_d(i0, i1, dim, dim);
637 FieldContainer<double> mb_x_x_d_d(i0, i1, dim, dim);
638 FieldContainer<double> mc_x_x_d_d(i0, i1, dim, dim);
639 FieldContainer<double> va_x_x_d(i0, i1, dim);
640 FieldContainer<double> vb_x_x_d(i0, i1, dim);
641 FieldContainer<double> vc_x_x_d(i0, i1, dim);
642 FieldContainer<double> vdot_x_x(i0, i1);
643 FieldContainer<double> vnorms_x_x(i0, i1);
644 FieldContainer<double> vnorms_x(i0);
645 double zero = INTREPID_TOL*100.0;
646
647 // fill with random numbers
648 for (int i=0; i<ma_x_x_d_d.size(); i++) {
649 ma_x_x_d_d[i] = Teuchos::ScalarTraits<double>::random();
650 }
651 for (int i=0; i<va_x_x_d.size(); i++) {
652 va_x_x_d[i] = Teuchos::ScalarTraits<double>::random();
653 }
654
655 *outStream << "\n************ Checking vectorNorm ************\n";
656
657 rst::vectorNorm(vnorms_x_x, va_x_x_d, NORM_TWO);
658 *outStream << va_x_x_d;
659 *outStream << vnorms_x_x;
660 if ( std::abs(rst::vectorNorm(&vnorms_x_x[0], vnorms_x_x.size(), NORM_TWO) -
661 rst::vectorNorm(&va_x_x_d[0], va_x_x_d.size(), NORM_TWO)) > zero) {
662 *outStream << "\n\nINCORRECT vectorNorm NORM_TWO\n\n";
663 errorFlag = -1000;
664 }
665
666 rst::vectorNorm(vnorms_x_x, va_x_x_d, NORM_ONE);
667 *outStream << va_x_x_d;
668 *outStream << vnorms_x_x;
669 if ( std::abs(rst::vectorNorm(&vnorms_x_x[0], vnorms_x_x.size(), NORM_ONE) -
670 rst::vectorNorm(&va_x_x_d[0], va_x_x_d.size(), NORM_ONE)) > zero) {
671 *outStream << "\n\nINCORRECT vectorNorm NORM_ONE\n\n";
672 errorFlag = -1000;
673 }
674
675 rst::vectorNorm(vnorms_x_x, va_x_x_d, NORM_INF);
676 *outStream << va_x_x_d;
677 *outStream << vnorms_x_x;
678 if ( std::abs(rst::vectorNorm(&vnorms_x_x[0], vnorms_x_x.size(), NORM_INF) -
679 rst::vectorNorm(&va_x_x_d[0], va_x_x_d.size(), NORM_INF)) > zero) {
680 *outStream << "\n\nINCORRECT vectorNorm NORM_INF\n\n";
681 errorFlag = -1000;
682 }
683
684 /******************************************/
685
686
687 *outStream << "\n************ Checking inverse, subtract, and vectorNorm ************\n";
688
689 rst::inverse(mb_x_x_d_d, ma_x_x_d_d); // B = inv(A)
690 rst::inverse(mc_x_x_d_d, mb_x_x_d_d); // C = inv(B) ~= A
691 *outStream << ma_x_x_d_d << mb_x_x_d_d << mc_x_x_d_d;
692
693 rst::subtract(&mc_x_x_d_d[0], &ma_x_x_d_d[0], ma_x_x_d_d.size()); // C = C - A ~= 0
694
695 if (rst::vectorNorm(&mc_x_x_d_d[0], mc_x_x_d_d.size(), NORM_ONE) > zero) {
696 *outStream << "\n\nINCORRECT inverse OR subtract OR vectorNorm\n\n";
697 errorFlag = -1000;
698 }
699
700 /******************************************/
701
702
703 *outStream << "\n********** Checking determinant **********\n";
704
705 FieldContainer<double> detA_x_x(i0, i1);
706 FieldContainer<double> detB_x_x(i0, i1);
707
708 rst::det(detA_x_x, ma_x_x_d_d);
709 rst::det(detB_x_x, mb_x_x_d_d);
710 *outStream << detA_x_x << detB_x_x;
711
712 if ( (rst::dot(&detA_x_x[0], &detB_x_x[0], detA_x_x.size()) - (double)(i0*i1)) > zero) {
713 *outStream << "\n\nINCORRECT det\n\n" ;
714 errorFlag = -1000;
715 }
716
717 *outStream << "\n det(A)*det(inv(A)) = " <<
718 rst::det(&ma_x_x_d_d[0], ma_x_x_d_d.dimension(3))*rst::det(&mb_x_x_d_d[0], mb_x_x_d_d.dimension(3))
719 << "\n";
720
721 if ( (rst::det(&ma_x_x_d_d[0], ma_x_x_d_d.dimension(3))*
722 rst::det(&mb_x_x_d_d[0], mb_x_x_d_d.dimension(3)) - (double)1) > zero) {
723 *outStream << "\n\nINCORRECT det\n\n" ;
724 errorFlag = -1000;
725 }
726
727 /******************************************/
728
729
730 *outStream << "\n************ Checking transpose and subtract ************\n";
731
732 rst::transpose(mb_x_x_d_d, ma_x_x_d_d); // B = A^T
733 rst::transpose(mc_x_x_d_d, mb_x_x_d_d); // C = B^T = A
734 *outStream << ma_x_x_d_d << mb_x_x_d_d << mc_x_x_d_d;
735
736 rst::subtract(&mc_x_x_d_d[0], &ma_x_x_d_d[0], ma_x_x_d_d.size()); // C = C - A = 0
737
738 if (rst::vectorNorm(&mc_x_x_d_d[0], mc_x_x_d_d.size(), NORM_ONE) > zero) {
739 *outStream << "\n\nINCORRECT transpose OR subtract OR vectorNorm\n\n" ;
740 errorFlag = -1000;
741 }
742
743 /******************************************/
744
745
746 *outStream << "\n************ Checking matvec, vectorNorm, subtract, and inverse ************\n";
747
748 rst::inverse(mb_x_x_d_d, ma_x_x_d_d); // B = inv(A)
749 rst::inverse(mc_x_x_d_d, mb_x_x_d_d); // C = inv(B) ~= A
750 rst::matvec(vb_x_x_d, ma_x_x_d_d, va_x_x_d); // b = A*a
751 rst::matvec(vc_x_x_d, mb_x_x_d_d, vb_x_x_d); // c = inv(A)*(A*a) ~= a
752 rst::subtract(vc_x_x_d, va_x_x_d); // c = c - a ~= 0
753 *outStream << vc_x_x_d;
754
755 rst::vectorNorm(vnorms_x_x, vc_x_x_d, NORM_ONE);
756 rst::vectorNorm(vnorms_x, vnorms_x_x, NORM_INF);
757 if (rst::vectorNorm(vnorms_x, NORM_TWO) > zero) {
758 *outStream << "\n\nINCORRECT matvec OR inverse OR subtract OR vectorNorm\n\n";
759 errorFlag = -1000;
760 }
761
762 /******************************************/
763
764
765 *outStream << "\n************ Checking add, subtract, absval, and scale ************\n";
766
767 double x = 1.234;
768 rst::add(vc_x_x_d, va_x_x_d, vb_x_x_d); // c = a + b
769 rst::subtract(vc_x_x_d, vb_x_x_d); // c = c - b = a
770 rst::scale(vb_x_x_d, vc_x_x_d, x); // b = c*x;
771 rst::scale(vc_x_x_d, vb_x_x_d, (1.0/x)); // c = b*(1/x) = a;
772 rst::subtract(vb_x_x_d, vc_x_x_d, va_x_x_d); // b = c - a ~= 0
773 rst::absval(vc_x_x_d, vb_x_x_d); // c = |b|
774 rst::scale(vb_x_x_d, vc_x_x_d, -1.0); // b = -c
775 rst::absval(vc_x_x_d, vb_x_x_d); // c = |b|
776 rst::add(vc_x_x_d, vb_x_x_d); // c = c + b === 0
777 *outStream << vc_x_x_d;
778
779 rst::vectorNorm(vnorms_x_x, vc_x_x_d, NORM_ONE);
780 rst::vectorNorm(vnorms_x, vnorms_x_x, NORM_INF);
781 if (rst::vectorNorm(vnorms_x, NORM_TWO) > (double)0) {
782 *outStream << "\n\nSign flips combined with std::abs might not be invertible on this platform!\n"
783 << "Potential IEEE compliance issues!\n\n";
784 if (rst::vectorNorm(vnorms_x, NORM_TWO) > zero) {
785 *outStream << "\n\nINCORRECT add OR subtract OR scale OR absval OR vectorNorm\n\n";
786 errorFlag = -1000;
787 }
788 }
789
790 /******************************************/
791
792
793 *outStream << "\n************ Checking dot and vectorNorm ************\n";
794
795 for (int i=0; i<va_x_x_d.size(); i++) {
796 va_x_x_d[i] = 2.0;
797 }
798
799 rst::dot(vdot_x_x, va_x_x_d, va_x_x_d); // dot = a'*a
800 *outStream << vdot_x_x;
801
802 rst::vectorNorm(vnorms_x, vdot_x_x, NORM_ONE);
803 if (rst::vectorNorm(vnorms_x, NORM_ONE) - (double)(4.0*dim*i0*i1) > zero) {
804 *outStream << "\n\nINCORRECT dot OR vectorNorm\n\n";
805 errorFlag = -1000;
806 }
807
808 /******************************************/
809
810 *outStream << "\n";
811 }
812
813 // one-dimensional base containers
814 for (int dim=3; dim>0; dim--) {
815 int i0=7;
816 FieldContainer<double> ma_x_d_d(i0, dim, dim);
817 FieldContainer<double> mb_x_d_d(i0, dim, dim);
818 FieldContainer<double> mc_x_d_d(i0, dim, dim);
819 FieldContainer<double> va_x_d(i0, dim);
820 FieldContainer<double> vb_x_d(i0, dim);
821 FieldContainer<double> vc_x_d(i0, dim);
822 FieldContainer<double> vdot_x(i0);
823 FieldContainer<double> vnorms_x(i0);
824 double zero = INTREPID_TOL*100.0;
825
826 // fill with random numbers
827 for (int i=0; i<ma_x_d_d.size(); i++) {
828 ma_x_d_d[i] = Teuchos::ScalarTraits<double>::random();
829 }
830 for (int i=0; i<va_x_d.size(); i++) {
831 va_x_d[i] = Teuchos::ScalarTraits<double>::random();
832 }
833
834 *outStream << "\n************ Checking vectorNorm ************\n";
835
836 rst::vectorNorm(vnorms_x, va_x_d, NORM_TWO);
837 *outStream << va_x_d;
838 *outStream << vnorms_x;
839 if ( std::abs(rst::vectorNorm(&vnorms_x[0], vnorms_x.size(), NORM_TWO) -
840 rst::vectorNorm(&va_x_d[0], va_x_d.size(), NORM_TWO)) > zero) {
841 *outStream << "\n\nINCORRECT vectorNorm NORM_TWO\n\n";
842 errorFlag = -1000;
843 }
844
845 rst::vectorNorm(vnorms_x, va_x_d, NORM_ONE);
846 *outStream << va_x_d;
847 *outStream << vnorms_x;
848 if ( std::abs(rst::vectorNorm(&vnorms_x[0], vnorms_x.size(), NORM_ONE) -
849 rst::vectorNorm(&va_x_d[0], va_x_d.size(), NORM_ONE)) > zero) {
850 *outStream << "\n\nINCORRECT vectorNorm NORM_ONE\n\n";
851 errorFlag = -1000;
852 }
853
854 rst::vectorNorm(vnorms_x, va_x_d, NORM_INF);
855 *outStream << va_x_d;
856 *outStream << vnorms_x;
857 if ( std::abs(rst::vectorNorm(&vnorms_x[0], vnorms_x.size(), NORM_INF) -
858 rst::vectorNorm(&va_x_d[0], va_x_d.size(), NORM_INF)) > zero) {
859 *outStream << "\n\nINCORRECT vectorNorm NORM_INF\n\n";
860 errorFlag = -1000;
861 }
862
863 /******************************************/
864
865
866 *outStream << "\n************ Checking inverse, subtract, and vectorNorm ************\n";
867
868 rst::inverse(mb_x_d_d, ma_x_d_d); // B = inv(A)
869 rst::inverse(mc_x_d_d, mb_x_d_d); // C = inv(B) ~= A
870 *outStream << ma_x_d_d << mb_x_d_d << mc_x_d_d;
871
872 rst::subtract(&mc_x_d_d[0], &ma_x_d_d[0], ma_x_d_d.size()); // C = C - A ~= 0
873
874 if (rst::vectorNorm(&mc_x_d_d[0], mc_x_d_d.size(), NORM_ONE) > zero) {
875 *outStream << "\n\nINCORRECT inverse OR subtract OR vectorNorm\n\n";
876 errorFlag = -1000;
877 }
878
879 /******************************************/
880
881
882 *outStream << "\n********** Checking determinant **********\n";
883
884 FieldContainer<double> detA_x(i0);
885 FieldContainer<double> detB_x(i0);
886
887 rst::det(detA_x, ma_x_d_d);
888 rst::det(detB_x, mb_x_d_d);
889 *outStream << detA_x << detB_x;
890
891 if ( (rst::dot(&detA_x[0], &detB_x[0], detA_x.size()) - (double)i0) > zero) {
892 *outStream << "\n\nINCORRECT det\n\n" ;
893 errorFlag = -1000;
894 }
895
896 *outStream << "\n det(A)*det(inv(A)) = " <<
897 rst::det(&ma_x_d_d[0], ma_x_d_d.dimension(2))*rst::det(&mb_x_d_d[0], mb_x_d_d.dimension(2))
898 << "\n";
899
900 if ( (rst::det(&ma_x_d_d[0], ma_x_d_d.dimension(2))*
901 rst::det(&mb_x_d_d[0], mb_x_d_d.dimension(2)) - (double)1) > zero) {
902 *outStream << "\n\nINCORRECT det\n\n" ;
903 errorFlag = -1000;
904 }
905
906 /******************************************/
907
908
909 *outStream << "\n************ Checking transpose and subtract ************\n";
910
911 rst::transpose(mb_x_d_d, ma_x_d_d); // B = A^T
912 rst::transpose(mc_x_d_d, mb_x_d_d); // C = B^T = A
913 *outStream << ma_x_d_d << mb_x_d_d << mc_x_d_d;
914
915 rst::subtract(&mc_x_d_d[0], &ma_x_d_d[0], ma_x_d_d.size()); // C = C - A = 0
916
917 if (rst::vectorNorm(&mc_x_d_d[0], mc_x_d_d.size(), NORM_ONE) > zero) {
918 *outStream << "\n\nINCORRECT transpose OR subtract OR vectorNorm\n\n" ;
919 errorFlag = -1000;
920 }
921
922 /******************************************/
923
924
925 *outStream << "\n************ Checking matvec, vectorNorm, subtract, and inverse ************\n";
926
927 rst::inverse(mb_x_d_d, ma_x_d_d); // B = inv(A)
928 rst::inverse(mc_x_d_d, mb_x_d_d); // C = inv(B) ~= A
929 rst::matvec(vb_x_d, ma_x_d_d, va_x_d); // b = A*a
930 rst::matvec(vc_x_d, mb_x_d_d, vb_x_d); // c = inv(A)*(A*a) ~= a
931 rst::subtract(vc_x_d, va_x_d); // c = c - a ~= 0
932 *outStream << vc_x_d;
933
934 rst::vectorNorm(vnorms_x, vc_x_d, NORM_ONE);
935 if (rst::vectorNorm(vnorms_x, NORM_TWO) > zero) {
936 *outStream << "\n\nINCORRECT matvec OR inverse OR subtract OR vectorNorm\n\n";
937 errorFlag = -1000;
938 }
939
940 /******************************************/
941
942
943 *outStream << "\n************ Checking add, subtract, absval, and scale ************\n";
944
945 double x = 1.234;
946 rst::add(vc_x_d, va_x_d, vb_x_d); // c = a + b
947 rst::subtract(vc_x_d, vb_x_d); // c = c - b = a
948 rst::scale(vb_x_d, vc_x_d, x); // b = c*x;
949 rst::scale(vc_x_d, vb_x_d, (1.0/x)); // c = b*(1/x) = a;
950 rst::subtract(vb_x_d, vc_x_d, va_x_d); // b = c - a ~= 0
951 rst::absval(vc_x_d, vb_x_d); // c = |b|
952 rst::scale(vb_x_d, vc_x_d, -1.0); // b = -c
953 rst::absval(vc_x_d, vb_x_d); // c = |b|
954 rst::add(vc_x_d, vb_x_d); // c = c + b === 0
955 *outStream << vc_x_d;
956
957 rst::vectorNorm(vnorms_x, vc_x_d, NORM_ONE);
958 if (rst::vectorNorm(vnorms_x, NORM_TWO) > (double)0) {
959 *outStream << "\n\nSign flips combined with std::abs might not be invertible on this platform!\n"
960 << "Potential IEEE compliance issues!\n\n";
961 if (rst::vectorNorm(vnorms_x, NORM_TWO) > zero) {
962 *outStream << "\n\nINCORRECT add OR subtract OR scale OR absval OR vectorNorm\n\n";
963 errorFlag = -1000;
964 }
965 }
966
967 /******************************************/
968
969
970 *outStream << "\n************ Checking dot and vectorNorm ************\n";
971
972 for (int i=0; i<va_x_d.size(); i++) {
973 va_x_d[i] = 2.0;
974 }
975 rst::dot(vdot_x, va_x_d, va_x_d); // dot = a'*a
976 *outStream << vdot_x;
977
978 if (rst::vectorNorm(vdot_x, NORM_ONE) - (double)(4.0*dim*i0) > zero) {
979 *outStream << "\n\nINCORRECT dot OR vectorNorm\n\n";
980 errorFlag = -1000;
981 }
982
983 /******************************************/
984
985 *outStream << "\n";
986 }
987 }
988 catch (const std::logic_error & err) {
989 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
990 *outStream << err.what() << '\n';
991 *outStream << "-------------------------------------------------------------------------------" << "\n\n";
992 errorFlag = -1000;
993 };
994
995 if (errorFlag != 0)
996 std::cout << "End Result: TEST FAILED\n";
997 else
998 std::cout << "End Result: TEST PASSED\n";
999
1000 // reset format state of std::cout
1001 std::cout.copyfmt(oldFormatState);
1002
1003 return errorFlag;
1004}
main
int main(int argc, char *argv[])
Performs a code-code comparison to FIAT for Nedelec bases on tets (values and curls)
Definition: test_01.cpp:65
Intrepid_FieldContainer.hpp
Header file for utility class to provide multidimensional containers.
Intrepid_RealSpaceTools.hpp
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D.
Intrepid::RealSpaceTools
Implementation of basic linear algebra functionality in Euclidean space.
Definition: Intrepid_RealSpaceTools.hpp:65
Generated by doxygen 1.9.4