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 
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 }
Implementation of basic linear algebra functionality in Euclidean space.
Header file for utility class to provide multidimensional containers.
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D...