Intrepid
test_01.cpp
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 
57 using namespace std;
58 using namespace Intrepid;
59 
60 int 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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 (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 }
Intrepid::RealSpaceTools
Implementation of basic linear algebra functionality in Euclidean space.
Definition: Intrepid_RealSpaceTools.hpp:68
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...
Generated by   doxygen 1.8.5