Intrepid
test_01_kokkos.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 #include <Kokkos_Core.hpp>
56 
57 
58 using namespace std;
59 using namespace Intrepid;
60 
61 int main(int argc, char *argv[]) {
62 
63  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
64 Kokkos::initialize();
65 
66  // This little trick lets us print to std::cout only if
67  // a (dummy) command-line argument is provided.
68  int iprint = argc - 1;
69  Teuchos::RCP<std::ostream> outStream;
70  Teuchos::oblackholestream bhs; // outputs nothing
71  if (iprint > 0)
72  outStream = Teuchos::rcp(&std::cout, false);
73  else
74  outStream = Teuchos::rcp(&bhs, false);
75 
76  // Save the format state of the original std::cout.
77  Teuchos::oblackholestream oldFormatState;
78  oldFormatState.copyfmt(std::cout);
79 
80  *outStream \
81  << "===============================================================================\n" \
82  << "| |\n" \
83  << "| Unit Test (RealSpaceTools) |\n" \
84  << "| |\n" \
85  << "| 1) Vector operations in 1D, 2D, and 3D real space |\n" \
86  << "| 2) Matrix / matrix-vector operations in 1D, 2D, and 3D real space |\n" \
87  << "| |\n" \
88  << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \
89  << "| Denis Ridzal (dridzal@sandia.gov). |\n" \
90  << "| |\n" \
91  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
92  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
93  << "| |\n" \
94  << "===============================================================================\n";
95 
96 
97  typedef RealSpaceTools<double> rst;
98 
99 
100  int errorFlag = 0;
101 #ifdef HAVE_INTREPID_DEBUG
102  int beginThrowNumber = Teuchos::TestForException_getThrowNumber();
103  int endThrowNumber = beginThrowNumber + 49;
104 #ifndef HAVE_INTREPID_DEBUG_INF_CHECK
105  endThrowNumber = beginThrowNumber + 44;
106 #endif
107 
108 #endif
109  *outStream \
110  << "\n"
111  << "===============================================================================\n"\
112  << "| TEST 1: vector exceptions |\n"\
113  << "===============================================================================\n";
114 
115  try{
116 
117  Kokkos::View<double**> a_2_2("a_2_2",2, 2);
118  Kokkos::View<double**> a_10_2("a_10_2",10, 2);
119  Kokkos::View<double**> a_10_3("a_10_3",10, 3);
120  Kokkos::View<double***> a_10_2_2("a_10_2_2",10, 2, 2);
121  Kokkos::View<double***> a_10_2_3("a_10_2_3",10, 2, 3);
122  Kokkos::View<double****> a_10_2_2_3("a_10_2_2_3",10, 2, 2, 3);
123 
124 #ifdef HAVE_INTREPID_DEBUG
125 
126  *outStream << "-> vector norm with multidimensional arrays:\n";
127 
128  try {
129  rst::vectorNorm(a_2_2, NORM_TWO);
130  }
131  catch (std::logic_error err) {
132  *outStream << "Expected Error ----------------------------------------------------------------\n";
133  *outStream << err.what() << '\n';
134  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
135  };
136  try {
137  rst::vectorNorm(a_10_2_2, a_10_2_2, NORM_TWO);
138  }
139  catch (std::logic_error err) {
140  *outStream << "Expected Error ----------------------------------------------------------------\n";
141  *outStream << err.what() << '\n';
142  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
143  };
144  try {
145  rst::vectorNorm(a_10_2_2, a_10_2_2_3, NORM_TWO);
146  }
147  catch (std::logic_error err) {
148  *outStream << "Expected Error ----------------------------------------------------------------\n";
149  *outStream << err.what() << '\n';
150  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
151  };
152  try {
153  rst::vectorNorm(a_10_3, a_10_2_2, NORM_TWO);
154  }
155  catch (std::logic_error err) {
156  *outStream << "Expected Error ----------------------------------------------------------------\n";
157  *outStream << err.what() << '\n';
158  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
159  };
160 
161  *outStream << "-> add with multidimensional arrays:\n";
162 
163  try {
164  rst::add(a_10_2_2, a_10_2, a_2_2);
165  }
166  catch (std::logic_error err) {
167  *outStream << "Expected Error ----------------------------------------------------------------\n";
168  *outStream << err.what() << '\n';
169  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
170  };
171  try {
172  rst::add(a_10_2_3, a_10_2_2, a_10_2_2);
173  }
174  catch (std::logic_error err) {
175  *outStream << "Expected Error ----------------------------------------------------------------\n";
176  *outStream << err.what() << '\n';
177  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
178  };
179  try {
180  rst::add(a_10_2_2, a_10_2_2_3);
181  }
182  catch (std::logic_error err) {
183  *outStream << "Expected Error ----------------------------------------------------------------\n";
184  *outStream << err.what() << '\n';
185  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
186  };
187  try {
188  rst::add(a_10_2_3, a_10_2_2);
189  }
190  catch (std::logic_error err) {
191  *outStream << "Expected Error ----------------------------------------------------------------\n";
192  *outStream << err.what() << '\n';
193  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
194  };
195 
196  *outStream << "-> subtract with multidimensional arrays:\n";
197 
198  try {
199  rst::subtract(a_10_2_2, a_10_2, a_2_2);
200  }
201  catch (std::logic_error err) {
202  *outStream << "Expected Error ----------------------------------------------------------------\n";
203  *outStream << err.what() << '\n';
204  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
205  };
206  try {
207  rst::subtract(a_10_2_3, a_10_2_2, a_10_2_2);
208  }
209  catch (std::logic_error err) {
210  *outStream << "Expected Error ----------------------------------------------------------------\n";
211  *outStream << err.what() << '\n';
212  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
213  };
214  try {
215  rst::subtract(a_10_2_2, a_10_2_2_3);
216  }
217  catch (std::logic_error err) {
218  *outStream << "Expected Error ----------------------------------------------------------------\n";
219  *outStream << err.what() << '\n';
220  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
221  };
222  try {
223  rst::subtract(a_10_2_3, a_10_2_2);
224  }
225  catch (std::logic_error err) {
226  *outStream << "Expected Error ----------------------------------------------------------------\n";
227  *outStream << err.what() << '\n';
228  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
229  };
230 
231  *outStream << "-> dot product norm with multidimensional arrays:\n";
232 
233  try {
234  rst::dot(a_10_2, a_10_2_2_3, a_10_2_2_3);
235  }
236  catch (std::logic_error err) {
237  *outStream << "Expected Error ----------------------------------------------------------------\n";
238  *outStream << err.what() << '\n';
239  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
240  };
241  try {
242  rst::dot(a_10_2, a_10_2_2, a_10_2_2_3);
243  }
244  catch (std::logic_error err) {
245  *outStream << "Expected Error ----------------------------------------------------------------\n";
246  *outStream << err.what() << '\n';
247  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
248  };
249  try {
250  rst::dot(a_10_2_2, a_10_2_2_3, a_10_2_2_3);
251  }
252  catch (std::logic_error err) {
253  *outStream << "Expected Error ----------------------------------------------------------------\n";
254  *outStream << err.what() << '\n';
255  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
256  };
257  try {
258  rst::dot(a_10_2, a_10_2_2, a_10_2_3);
259  }
260  catch (std::logic_error err) {
261  *outStream << "Expected Error ----------------------------------------------------------------\n";
262  *outStream << err.what() << '\n';
263  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
264  };
265  try {
266  rst::dot(a_10_3, a_10_2_3, a_10_2_3);
267  }
268  catch (std::logic_error err) {
269  *outStream << "Expected Error ----------------------------------------------------------------\n";
270  *outStream << err.what() << '\n';
271  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
272  };
273 
274  *outStream << "-> absolute value with multidimensional arrays:\n";
275 
276  try {
277  rst::absval(a_10_3, a_10_2_3);
278  }
279  catch (std::logic_error err) {
280  *outStream << "Expected Error ----------------------------------------------------------------\n";
281  *outStream << err.what() << '\n';
282  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
283  };
284  try {
285  rst::absval(a_10_2_2, a_10_2_3);
286  }
287  catch (std::logic_error err) {
288  *outStream << "Expected Error ----------------------------------------------------------------\n";
289  *outStream << err.what() << '\n';
290  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
291  };
292 #endif
293 
294  }
295  catch (std::logic_error err) {
296  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
297  *outStream << err.what() << '\n';
298  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
299  errorFlag = -1000;
300  };
301 
302 
303  *outStream \
304  << "\n"
305  << "===============================================================================\n"\
306  << "| TEST 2: matrix / matrix-vector exceptions |\n"\
307  << "===============================================================================\n";
308 
309  try{
310 
311  Kokkos::View<double*****> a_10_1_2_3_4("a_10_1_2_3_4",10, 1, 2, 3, 4);
312  Kokkos::View<double*****> b_10_1_2_3_4("b_10_1_2_3_4",10, 1, 2, 3, 4);
313  Kokkos::View<double*> a_10("a_10",10);
314  Kokkos::View<double*> a_9("a_9",9);
315  Kokkos::View<double*> b_10("b_10",10);
316  Kokkos::View<double****> a_10_15_4_4("a_10_15_4_4",10, 15, 4, 4);
317  Kokkos::View<double****> b_10_15_4_4("b_10_15_4_4",10, 15, 4, 4);
318  Kokkos::View<double***> a_10_2_2("a_10_2_2",10, 2, 2);
319  Kokkos::View<double***> a_10_2_3("a_10_2_3",10, 2, 3);
320  Kokkos::View<double***> b_10_2_3("b_10_2_3",10, 2, 3);
321 
322  Kokkos::View<double**> a_1_1("a_1_1",1, 1);
323  Kokkos::View<double**> b_1_1("b_1_1",1, 1);
324  Kokkos::View<double**> a_2_2("a_2_2",2, 2);
325  Kokkos::View<double**> b_2_2("b_2_2",2, 2);
326  Kokkos::View<double**> a_3_3("a_3_3",3, 3);
327  Kokkos::View<double**> b_3_3("b_3_3",3, 3);
328  Kokkos::View<double**> a_2_3("a_2_3",2, 3);
329  Kokkos::View<double**> a_4_4("a_4_4",4, 4);
330 
331  Kokkos::View<double****> a_10_15_1_1("a_10_15_1_1",10, 15, 1, 1);
332  Kokkos::View<double****> b_10_15_1_1("b_10_15_1_1",10, 15, 1, 1);
333  Kokkos::View<double****> a_10_15_2_2("a_10_15_2_2",10, 15, 2, 2);
334  Kokkos::View<double****> b_10_15_2_2("b_10_15_2_2",10, 15, 2, 2);
335  Kokkos::View<double****> a_10_15_3_3("a_10_15_3_3",10, 15, 3, 3);
336  Kokkos::View<double****> a_10_15_3_2("a_10_15_3_2",10, 15, 3, 2);
337  Kokkos::View<double****> b_10_15_3_3("b_10_15_3_3",10, 15, 3, 3);
338  Kokkos::View<double**> b_10_14("b_10_14",10, 14);
339  Kokkos::View<double**> b_10_15("b_10_15",10, 15);
340  Kokkos::View<double***> b_10_14_3("b_10_14_3",10, 14, 3);
341  Kokkos::View<double***> b_10_15_3("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,0) = 1.0;
421  a_3_3(0,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  *outStream \
623  << "\n"
624  << "===============================================================================\n"\
625  << "| TEST 2: correctness of math operations |\n"\
626  << "===============================================================================\n";
627 
628  outStream->precision(20);
629 
630  try{
631 
632  // two-dimensional base containers
633  for (int dim=3; dim>0; dim--) {
634  int i0=4, i1=5;
635  Kokkos::View<double****> ma_x_x_d_d("ma_x_x_d_d",i0, i1, dim, dim);
636  Kokkos::View<double****> mb_x_x_d_d("mb_x_x_d_d",i0, i1, dim, dim);
637  Kokkos::View<double****> mc_x_x_d_d("mc_x_x_d_d",i0, i1, dim, dim);
638  Kokkos::View<double***> va_x_x_d("va_x_x_d",i0, i1, dim);
639  Kokkos::View<double***> vb_x_x_d("vb_x_x_d",i0, i1, dim);
640  Kokkos::View<double***> vc_x_x_d("vc_x_x_d",i0, i1, dim);
641  Kokkos::View<double**> vdot_x_x("vdot_x_x",i0, i1);
642  Kokkos::View<double**> vnorms_x_x("vnorms_x_x",i0, i1);
643  Kokkos::View<double*> vnorms_x("vnorms_x",i0);
644  double zero = INTREPID_TOL*100.0;
645 
646  // fill with random numbers
647  for (unsigned int i=0; i<ma_x_x_d_d.dimension(0); i++) {
648  for (unsigned int j=0; j<ma_x_x_d_d.dimension(1); j++) {
649  for (unsigned int k=0; k<ma_x_x_d_d.dimension(2); k++) {
650  for (unsigned int l=0; l<ma_x_x_d_d.dimension(3); l++) {
651  ma_x_x_d_d(i,j,k,l) = Teuchos::ScalarTraits<double>::random();
652 
653  }
654  }
655  }
656  }
657  for (unsigned int i=0; i<va_x_x_d.dimension(0); i++) {
658  for (unsigned int j=0; j<va_x_x_d.dimension(1); j++) {
659  for (unsigned int k=0; k<va_x_x_d.dimension(2); k++) {
660  va_x_x_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
661  }
662  }
663  }
664 
665  *outStream << "\n************ Checking vectorNorm ************\n";
666 
667  rst::vectorNorm(vnorms_x_x, va_x_x_d, NORM_TWO);
668  *outStream << "va_x_x_d";
669  *outStream << "vnorms_x_x";
670  if ( std::abs(rst::vectorNorm(vnorms_x_x, NORM_TWO) -
671  rst::vectorNorm(va_x_x_d, NORM_TWO)) > zero) {
672  *outStream << "\n\nINCORRECT vectorNorm NORM_TWO\n\n";
673  errorFlag = -1000;
674  }
675 
676  rst::vectorNorm(vnorms_x_x, va_x_x_d, NORM_ONE);
677  *outStream << "va_x_x_d";
678  *outStream << "vnorms_x_x";
679  if ( std::abs(rst::vectorNorm(vnorms_x_x, NORM_ONE) -
680  rst::vectorNorm(va_x_x_d, NORM_ONE)) > zero) {
681  *outStream << "\n\nINCORRECT vectorNorm NORM_ONE\n\n";
682  errorFlag = -1000;
683  }
684  rst::vectorNorm(vnorms_x_x, va_x_x_d, NORM_INF);
685  *outStream << "va_x_x_d";
686  *outStream << "vnorms_x_x";
687  if ( std::abs(rst::vectorNorm(vnorms_x_x, NORM_INF) -
688  rst::vectorNorm(va_x_x_d, NORM_INF)) > zero) {
689  *outStream << "\n\nINCORRECT vectorNorm NORM_INF\n\n";
690  errorFlag = -1000;
691  }
692 
693  /******************************************/
694 
695 
696  *outStream << "\n************ Checking inverse, subtract, and vectorNorm ************\n";
697 
698  rst::inverse(mb_x_x_d_d, ma_x_x_d_d); // B = inv(A)
699  rst::inverse(mc_x_x_d_d, mb_x_x_d_d); // C = inv(B) ~= A
700  *outStream << "ma_x_x_d_d" << "mb_x_x_d_d" << "mc_x_x_d_d";
701 
702  rst::subtract(mc_x_x_d_d, ma_x_x_d_d); // C = C - A ~= 0
703 
704  if (rst::vectorNorm(mc_x_x_d_d, NORM_ONE) > zero) {
705  *outStream << "\n\nINCORRECT inverse OR subtract OR vectorNorm\n\n";
706  errorFlag = -1000;
707  }
708 
709  /******************************************/
710 
711 
712  *outStream << "\n********** Checking determinant **********\n";
713 
714  Kokkos::View<double**> detA_x_x("detA_x_x",i0, i1);
715  Kokkos::View<double**> detB_x_x("detB_x_x",i0, i1);
716  Kokkos::View<double*> detA_x_xlinear("detA_x_xlinear",i0*i1);
717  Kokkos::View<double*> detB_x_xlinear("detB_x_xlinear",i0*i1);
718 
719  rst::det(detA_x_x, ma_x_x_d_d);
720  rst::det(detB_x_x, mb_x_x_d_d);
721  *outStream << "detA_x_x" << "detB_x_x";
722  for(int i=0;i<i0;i++){
723  for(int j=0;j<i1;j++){
724  detA_x_xlinear(i1*i+j)=detA_x_x(i,j);
725  detB_x_xlinear(i1*i+j)=detB_x_x(i,j);
726  }
727  }
728  if ( (rst::dot(detA_x_xlinear, detB_x_xlinear) - (double)(i0*i1)) > zero) {
729  *outStream << "\n\nINCORRECT det\n\n" ;
730  errorFlag = -1000;
731  }
732  /* Kokkos::View<double*> ma_x_x_d_dlinear("ma_x_x_d_d",i0*i1*dim*dim);
733  Kokkos::View<double*> mb_x_x_d_dlinear("mb_x_x_d_d",i0*i1*dim*dim);
734  for(int i=0;i<i0;i++){
735  for(int j=0;j<i1;j++){
736  for(int k=0;k<dim;k++){
737  for(int l=0;l<dim;l++){
738  ma_x_x_d_dlinear(i*i1*dim*dim+j*dim*dim+k*dim+l)=ma_x_x_d_d(i,j,k,l);
739  mb_x_x_d_dlinear(i*i1*dim*dim+j*dim*dim+k*dim+l)=mb_x_x_d_d(i,j,k,l);
740  }
741  }
742  }
743  }
744  *outStream << "\n det(A)*det(inv(A)) = " <<
745  rst::det(ma_x_x_d_dlinear)*rst::det(mb_x_x_d_dlinear)
746  << "\n";
747 
748 
749  if ( (rst::det(ma_x_x_d_dlinear)*
750  rst::det(mb_x_x_d_dlinear) - (double)1) > zero) {
751  *outStream << "\n\nINCORRECT det\n\n" ;
752  errorFlag = -1000;
753  }*/
754 
755  /******************************************/
756 
757 
758  *outStream << "\n************ Checking transpose and subtract ************\n";
759 
760  rst::transpose(mb_x_x_d_d, ma_x_x_d_d); // B = A^T
761  rst::transpose(mc_x_x_d_d, mb_x_x_d_d); // C = B^T = A
762  *outStream << "ma_x_x_d_d" << "mb_x_x_d_d" << "mc_x_x_d_d";
763 
764  rst::subtract(mc_x_x_d_d, ma_x_x_d_d); // C = C - A = 0
765 
766  if (rst::vectorNorm(mc_x_x_d_d, NORM_ONE) > zero) {
767  *outStream << "\n\nINCORRECT transpose OR subtract OR vectorNorm\n\n" ;
768  errorFlag = -1000;
769  }
770 
771  /******************************************/
772 
773 
774  *outStream << "\n************ Checking matvec, vectorNorm, subtract, and inverse ************\n";
775 
776  rst::inverse(mb_x_x_d_d, ma_x_x_d_d); // B = inv(A)
777  rst::inverse(mc_x_x_d_d, mb_x_x_d_d); // C = inv(B) ~= A
778  rst::matvec(vb_x_x_d, ma_x_x_d_d, va_x_x_d); // b = A*a
779  rst::matvec(vc_x_x_d, mb_x_x_d_d, vb_x_x_d); // c = inv(A)*(A*a) ~= a
780  rst::subtract(vc_x_x_d, va_x_x_d); // c = c - a ~= 0
781  *outStream << "vc_x_x_d";
782 
783  rst::vectorNorm(vnorms_x_x, vc_x_x_d, NORM_ONE);
784  rst::vectorNorm(vnorms_x, vnorms_x_x, NORM_INF);
785  if (rst::vectorNorm(vnorms_x, NORM_TWO) > zero) {
786  *outStream << "\n\nINCORRECT matvec OR inverse OR subtract OR vectorNorm\n\n";
787  errorFlag = -1000;
788  }
789 
790  /******************************************/
791 
792 
793  *outStream << "\n************ Checking add, subtract, absval, and scale ************\n";
794 
795  double x = 1.234;
796  rst::add(vc_x_x_d, va_x_x_d, vb_x_x_d); // c = a + b
797  rst::subtract(vc_x_x_d, vb_x_x_d); // c = c - b = a
798  rst::scale(vb_x_x_d, vc_x_x_d, x); // b = c*x;
799  rst::scale(vc_x_x_d, vb_x_x_d, (1.0/x)); // c = b*(1/x) = a;
800  rst::subtract(vb_x_x_d, vc_x_x_d, va_x_x_d); // b = c - a ~= 0
801  rst::absval(vc_x_x_d, vb_x_x_d); // c = |b|
802  rst::scale(vb_x_x_d, vc_x_x_d, -1.0); // b = -c
803  rst::absval(vc_x_x_d, vb_x_x_d); // c = |b|
804  rst::add(vc_x_x_d, vb_x_x_d); // c = c + b === 0
805  *outStream << "vc_x_x_d";
806 
807  rst::vectorNorm(vnorms_x_x, vc_x_x_d, NORM_ONE);
808  rst::vectorNorm(vnorms_x, vnorms_x_x, NORM_INF);
809  if (rst::vectorNorm(vnorms_x, NORM_TWO) > (double)0) {
810  *outStream << "\n\nSign flips combined with std::abs might not be invertible on this platform!\n"
811  << "Potential IEEE compliance issues!\n\n";
812  if (rst::vectorNorm(vnorms_x, NORM_TWO) > zero) {
813  *outStream << "\n\nINCORRECT add OR subtract OR scale OR absval OR vectorNorm\n\n";
814  errorFlag = -1000;
815  }
816  }
817 
818  /******************************************/
819 
820 
821  *outStream << "\n************ Checking dot and vectorNorm ************\n";
822  for (unsigned int i=0; i<va_x_x_d.dimension(0); i++) {
823  for (unsigned int j=0; j<va_x_x_d.dimension(1); j++) {
824  for (unsigned int k=0; k<va_x_x_d.dimension(2); k++) {
825  va_x_x_d(i,j,k) = 2.0;
826  }
827  }
828  }
829 
830 
831  rst::dot(vdot_x_x, va_x_x_d, va_x_x_d); // dot = a'*a
832  *outStream << "vdot_x_x";
833 
834  rst::vectorNorm(vnorms_x, vdot_x_x, NORM_ONE);
835  if (rst::vectorNorm(vnorms_x, NORM_ONE) - (double)(4.0*dim*i0*i1) > zero) {
836  *outStream << "\n\nINCORRECT dot OR vectorNorm\n\n";
837  errorFlag = -1000;
838  }
839 
840  /******************************************/
841 
842  *outStream << "\n";
843  }
844 
845  // one-dimensional base containers
846  for (int dim=3; dim>0; dim--) {
847  int i0=7;
848  Kokkos::View<double***> ma_x_d_d("ma_x_d_d",i0, dim, dim);
849  Kokkos::View<double***> mb_x_d_d("mb_x_d_d",i0, dim, dim);
850  Kokkos::View<double***> mc_x_d_d("mc_x_d_d",i0, dim, dim);
851  Kokkos::View<double**> va_x_d("va_x_d",i0, dim);
852  Kokkos::View<double**> vb_x_d("vb_x_d",i0, dim);
853  Kokkos::View<double**> vc_x_d("vc_x_d",i0, dim);
854  Kokkos::View<double*> vdot_x("vdot_x",i0);
855  Kokkos::View<double*> vnorms_x("vnorms_x",i0);
856  double zero = INTREPID_TOL*100.0;
857 
858  // fill with random numbers
859 
860  for (unsigned int i=0; i<ma_x_d_d.dimension(0); i++) {
861  for (unsigned int j=0; j<ma_x_d_d.dimension(1); j++) {
862  for (unsigned int k=0; k<ma_x_d_d.dimension(2); k++) {
863  ma_x_d_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
864 
865  }
866  }
867  }
868  for (unsigned int i=0; i<va_x_d.dimension(0); i++) {
869  for (unsigned int j=0; j<va_x_d.dimension(1); j++) {
870  va_x_d(i,j) = Teuchos::ScalarTraits<double>::random();
871 
872  }
873  }
874 
875 
876  *outStream << "\n************ Checking vectorNorm ************\n";
877 
878  rst::vectorNorm(vnorms_x, va_x_d, NORM_TWO);
879  *outStream << "va_x_d";
880  *outStream << "vnorms_x";
881  if ( std::abs(rst::vectorNorm(vnorms_x, NORM_TWO) -
882  rst::vectorNorm(va_x_d, NORM_TWO)) > zero) {
883  *outStream << "\n\nINCORRECT vectorNorm NORM_TWO\n\n";
884  errorFlag = -1000;
885  }
886 
887  rst::vectorNorm(vnorms_x, va_x_d, NORM_ONE);
888  *outStream << "va_x_d";
889  *outStream << "vnorms_x";
890  if ( std::abs(rst::vectorNorm(vnorms_x, NORM_ONE) -
891  rst::vectorNorm(va_x_d, NORM_ONE)) > zero) {
892  *outStream << "\n\nINCORRECT vectorNorm NORM_ONE\n\n";
893  errorFlag = -1000;
894  }
895 
896  rst::vectorNorm(vnorms_x, va_x_d, NORM_INF);
897  *outStream << "va_x_d";
898  *outStream << "vnorms_x";
899  if ( std::abs(rst::vectorNorm(vnorms_x, NORM_INF) -
900  rst::vectorNorm(va_x_d, NORM_INF)) > zero) {
901  *outStream << "\n\nINCORRECT vectorNorm NORM_INF\n\n";
902  errorFlag = -1000;
903  }
904 
905  /******************************************/
906 
907 
908  *outStream << "\n************ Checking inverse, subtract, and vectorNorm ************\n";
909 
910  rst::inverse(mb_x_d_d, ma_x_d_d); // B = inv(A)
911  rst::inverse(mc_x_d_d, mb_x_d_d); // C = inv(B) ~= A
912  *outStream << "ma_x_d_d" << "mb_x_d_d" << "mc_x_d_d";
913 
914  rst::subtract(mc_x_d_d, ma_x_d_d); // C = C - A ~= 0
915 
916  if (rst::vectorNorm(mc_x_d_d, NORM_ONE) > zero) {
917  *outStream << "\n\nINCORRECT inverse OR subtract OR vectorNorm\n\n";
918  errorFlag = -1000;
919  }
920 
921  /******************************************/
922 
923 
924  *outStream << "\n********** Checking determinant **********\n";
925 
926  FieldContainer<double> detA_x(i0);
927  FieldContainer<double> detB_x(i0);
928 
929  rst::det(detA_x, ma_x_d_d);
930  rst::det(detB_x, mb_x_d_d);
931  *outStream << "detA_x" << "detB_x";
932 
933  if ( (rst::dot(detA_x, detB_x) - (double)i0) > zero) {
934  *outStream << "\n\nINCORRECT det\n\n" ;
935  errorFlag = -1000;
936  }
937 /*
938  *outStream << "\n det(A)*det(inv(A)) = " <<
939  rst::det(ma_x_d_d)*rst::det(mb_x_d_d)
940  << "\n";
941 
942  if ( (rst::det(ma_x_d_d)*
943  rst::det(mb_x_d_d) - (double)1) > zero) {
944  *outStream << "\n\nINCORRECT det\n\n" ;
945  errorFlag = -1000;
946  }
947 */
948  /******************************************/
949 
950 
951  *outStream << "\n************ Checking transpose and subtract ************\n";
952 
953  rst::transpose(mb_x_d_d, ma_x_d_d); // B = A^T
954  rst::transpose(mc_x_d_d, mb_x_d_d); // C = B^T = A
955  *outStream << "ma_x_d_d" << "mb_x_d_d" << "mc_x_d_d";
956 
957  rst::subtract(mc_x_d_d, ma_x_d_d); // C = C - A = 0
958 
959  if (rst::vectorNorm(mc_x_d_d, NORM_ONE) > zero) {
960  *outStream << "\n\nINCORRECT transpose OR subtract OR vectorNorm\n\n" ;
961  errorFlag = -1000;
962  }
963 
964  /******************************************/
965 
966 
967  *outStream << "\n************ Checking matvec, vectorNorm, subtract, and inverse ************\n";
968 
969  rst::inverse(mb_x_d_d, ma_x_d_d); // B = inv(A)
970  rst::inverse(mc_x_d_d, mb_x_d_d); // C = inv(B) ~= A
971  rst::matvec(vb_x_d, ma_x_d_d, va_x_d); // b = A*a
972  rst::matvec(vc_x_d, mb_x_d_d, vb_x_d); // c = inv(A)*(A*a) ~= a
973  rst::subtract(vc_x_d, va_x_d); // c = c - a ~= 0
974  *outStream << "vc_x_d";
975 
976  rst::vectorNorm(vnorms_x, vc_x_d, NORM_ONE);
977  if (rst::vectorNorm(vnorms_x, NORM_TWO) > zero) {
978  *outStream << "\n\nINCORRECT matvec OR inverse OR subtract OR vectorNorm\n\n";
979  errorFlag = -1000;
980  }
981 
982  /******************************************/
983 
984 
985  *outStream << "\n************ Checking add, subtract, absval, and scale ************\n";
986 
987  double x = 1.234;
988  rst::add(vc_x_d, va_x_d, vb_x_d); // c = a + b
989  rst::subtract(vc_x_d, vb_x_d); // c = c - b = a
990  rst::scale(vb_x_d, vc_x_d, x); // b = c*x;
991  rst::scale(vc_x_d, vb_x_d, (1.0/x)); // c = b*(1/x) = a;
992  rst::subtract(vb_x_d, vc_x_d, va_x_d); // b = c - a ~= 0
993  rst::absval(vc_x_d, vb_x_d); // c = |b|
994  rst::scale(vb_x_d, vc_x_d, -1.0); // b = -c
995  rst::absval(vc_x_d, vb_x_d); // c = |b|
996  rst::add(vc_x_d, vb_x_d); // c = c + b === 0
997  *outStream << "vc_x_d";
998 
999  rst::vectorNorm(vnorms_x, vc_x_d, NORM_ONE);
1000  if (rst::vectorNorm(vnorms_x, NORM_TWO) > (double)0) {
1001  *outStream << "\n\nSign flips combined with std::abs might not be invertible on this platform!\n"
1002  << "Potential IEEE compliance issues!\n\n";
1003  if (rst::vectorNorm(vnorms_x, NORM_TWO) > zero) {
1004  *outStream << "\n\nINCORRECT add OR subtract OR scale OR absval OR vectorNorm\n\n";
1005  errorFlag = -1000;
1006  }
1007  }
1008 
1009  /******************************************/
1010 
1011 
1012  *outStream << "\n************ Checking dot and vectorNorm ************\n";
1013  for (unsigned int i=0; i<va_x_d.dimension(0); i++) {
1014  for (unsigned int j=0; j<va_x_d.dimension(1); j++) {
1015  va_x_d(i,j) = 2.0;
1016  }
1017  }
1018  rst::dot(vdot_x, va_x_d, va_x_d); // dot = a'*a
1019  *outStream << "vdot_x";
1020 
1021  if (rst::vectorNorm(vdot_x, NORM_ONE) - (double)(4.0*dim*i0) > zero) {
1022  *outStream << "\n\nINCORRECT dot OR vectorNorm\n\n";
1023  errorFlag = -1000;
1024  }
1025 
1026  /******************************************/
1027 
1028  *outStream << "\n";
1029  }
1030  }
1031  catch (std::logic_error err) {
1032  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
1033  *outStream << err.what() << '\n';
1034  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
1035  errorFlag = -1000;
1036  };
1037  if (errorFlag != 0)
1038  std::cout << "End Result: TEST FAILED\n";
1039  else
1040  std::cout << "End Result: TEST PASSED\n";
1041 
1042  // reset format state of std::cout
1043  std::cout.copyfmt(oldFormatState);
1044  Kokkos::finalize();
1045 
1046  return errorFlag;
1047 }
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