|
Intrepid
|
00001 // @HEADER 00002 // ************************************************************************ 00003 // 00004 // Intrepid Package 00005 // Copyright (2007) Sandia Corporation 00006 // 00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive 00008 // license for use of this work by or on behalf of the U.S. Government. 00009 // 00010 // Redistribution and use in source and binary forms, with or without 00011 // modification, are permitted provided that the following conditions are 00012 // met: 00013 // 00014 // 1. Redistributions of source code must retain the above copyright 00015 // notice, this list of conditions and the following disclaimer. 00016 // 00017 // 2. Redistributions in binary form must reproduce the above copyright 00018 // notice, this list of conditions and the following disclaimer in the 00019 // documentation and/or other materials provided with the distribution. 00020 // 00021 // 3. Neither the name of the Corporation nor the names of the 00022 // contributors may be used to endorse or promote products derived from 00023 // this software without specific prior written permission. 00024 // 00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY 00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE 00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00036 // 00037 // Questions? Contact Pavel Bochev (pbboche@sandia.gov) 00038 // Denis Ridzal (dridzal@sandia.gov), or 00039 // Kara Peterson (kjpeter@sandia.gov) 00040 // 00041 // ************************************************************************ 00042 // @HEADER 00043 00049 #include "Intrepid_RealSpaceTools.hpp" 00050 #include "Intrepid_FieldContainer.hpp" 00051 #include "Teuchos_oblackholestream.hpp" 00052 #include "Teuchos_RCP.hpp" 00053 #include "Teuchos_ScalarTraits.hpp" 00054 #include "Teuchos_GlobalMPISession.hpp" 00055 00056 00057 using namespace std; 00058 using namespace Intrepid; 00059 00060 int main(int argc, char *argv[]) { 00061 00062 Teuchos::GlobalMPISession mpiSession(&argc, &argv); 00063 00064 // This little trick lets us print to std::cout only if 00065 // a (dummy) command-line argument is provided. 00066 int iprint = argc - 1; 00067 Teuchos::RCP<std::ostream> outStream; 00068 Teuchos::oblackholestream bhs; // outputs nothing 00069 if (iprint > 0) 00070 outStream = Teuchos::rcp(&std::cout, false); 00071 else 00072 outStream = Teuchos::rcp(&bhs, false); 00073 00074 // Save the format state of the original std::cout. 00075 Teuchos::oblackholestream oldFormatState; 00076 oldFormatState.copyfmt(std::cout); 00077 00078 *outStream \ 00079 << "===============================================================================\n" \ 00080 << "| |\n" \ 00081 << "| Unit Test (RealSpaceTools) |\n" \ 00082 << "| |\n" \ 00083 << "| 1) Vector operations in 1D, 2D, and 3D real space |\n" \ 00084 << "| 2) Matrix / matrix-vector operations in 1D, 2D, and 3D real space |\n" \ 00085 << "| |\n" \ 00086 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \ 00087 << "| Denis Ridzal (dridzal@sandia.gov). |\n" \ 00088 << "| |\n" \ 00089 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \ 00090 << "| Trilinos website: http://trilinos.sandia.gov |\n" \ 00091 << "| |\n" \ 00092 << "===============================================================================\n"; 00093 00094 00095 typedef RealSpaceTools<double> rst; 00096 00097 00098 int errorFlag = 0; 00099 #ifdef HAVE_INTREPID_DEBUG 00100 int beginThrowNumber = Teuchos::TestForException_getThrowNumber(); 00101 int endThrowNumber = beginThrowNumber + 50; 00102 #ifndef HAVE_INTREPID_DEBUG_INF_CHECK 00103 endThrowNumber = beginThrowNumber + 45; 00104 #endif 00105 00106 #endif 00107 00108 *outStream \ 00109 << "\n" 00110 << "===============================================================================\n"\ 00111 << "| TEST 1: vector exceptions |\n"\ 00112 << "===============================================================================\n"; 00113 00114 try{ 00115 00116 FieldContainer<double> a_2_2(2, 2); 00117 FieldContainer<double> a_10_2(10, 2); 00118 FieldContainer<double> a_10_3(10, 3); 00119 FieldContainer<double> a_10_2_2(10, 2, 2); 00120 FieldContainer<double> a_10_2_3(10, 2, 3); 00121 FieldContainer<double> a_10_2_2_3(10, 2, 2, 3); 00122 00123 #ifdef HAVE_INTREPID_DEBUG 00124 00125 *outStream << "-> vector norm with multidimensional arrays:\n"; 00126 00127 try { 00128 rst::vectorNorm(a_2_2, NORM_TWO); 00129 } 00130 catch (std::logic_error err) { 00131 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00132 *outStream << err.what() << '\n'; 00133 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00134 }; 00135 try { 00136 rst::vectorNorm(a_10_2_2, a_10_2_2, NORM_TWO); 00137 } 00138 catch (std::logic_error err) { 00139 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00140 *outStream << err.what() << '\n'; 00141 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00142 }; 00143 try { 00144 rst::vectorNorm(a_10_2_2, a_10_2_2_3, NORM_TWO); 00145 } 00146 catch (std::logic_error err) { 00147 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00148 *outStream << err.what() << '\n'; 00149 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00150 }; 00151 try { 00152 rst::vectorNorm(a_10_3, a_10_2_2, NORM_TWO); 00153 } 00154 catch (std::logic_error err) { 00155 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00156 *outStream << err.what() << '\n'; 00157 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00158 }; 00159 00160 *outStream << "-> add with multidimensional arrays:\n"; 00161 00162 try { 00163 rst::add(a_10_2_2, a_10_2, a_2_2); 00164 } 00165 catch (std::logic_error err) { 00166 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00167 *outStream << err.what() << '\n'; 00168 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00169 }; 00170 try { 00171 rst::add(a_10_2_3, a_10_2_2, a_10_2_2); 00172 } 00173 catch (std::logic_error err) { 00174 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00175 *outStream << err.what() << '\n'; 00176 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00177 }; 00178 try { 00179 rst::add(a_10_2_2, a_10_2_2_3); 00180 } 00181 catch (std::logic_error err) { 00182 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00183 *outStream << err.what() << '\n'; 00184 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00185 }; 00186 try { 00187 rst::add(a_10_2_3, a_10_2_2); 00188 } 00189 catch (std::logic_error err) { 00190 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00191 *outStream << err.what() << '\n'; 00192 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00193 }; 00194 00195 *outStream << "-> subtract with multidimensional arrays:\n"; 00196 00197 try { 00198 rst::subtract(a_10_2_2, a_10_2, a_2_2); 00199 } 00200 catch (std::logic_error err) { 00201 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00202 *outStream << err.what() << '\n'; 00203 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00204 }; 00205 try { 00206 rst::subtract(a_10_2_3, a_10_2_2, a_10_2_2); 00207 } 00208 catch (std::logic_error err) { 00209 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00210 *outStream << err.what() << '\n'; 00211 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00212 }; 00213 try { 00214 rst::subtract(a_10_2_2, a_10_2_2_3); 00215 } 00216 catch (std::logic_error err) { 00217 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00218 *outStream << err.what() << '\n'; 00219 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00220 }; 00221 try { 00222 rst::subtract(a_10_2_3, a_10_2_2); 00223 } 00224 catch (std::logic_error err) { 00225 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00226 *outStream << err.what() << '\n'; 00227 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00228 }; 00229 00230 *outStream << "-> dot product norm with multidimensional arrays:\n"; 00231 00232 try { 00233 rst::dot(a_10_2, a_10_2_2_3, a_10_2_2_3); 00234 } 00235 catch (std::logic_error err) { 00236 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00237 *outStream << err.what() << '\n'; 00238 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00239 }; 00240 try { 00241 rst::dot(a_10_2, a_10_2_2, a_10_2_2_3); 00242 } 00243 catch (std::logic_error err) { 00244 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00245 *outStream << err.what() << '\n'; 00246 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00247 }; 00248 try { 00249 rst::dot(a_10_2_2, a_10_2_2_3, a_10_2_2_3); 00250 } 00251 catch (std::logic_error err) { 00252 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00253 *outStream << err.what() << '\n'; 00254 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00255 }; 00256 try { 00257 rst::dot(a_10_2, a_10_2_2, a_10_2_3); 00258 } 00259 catch (std::logic_error err) { 00260 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00261 *outStream << err.what() << '\n'; 00262 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00263 }; 00264 try { 00265 rst::dot(a_10_3, a_10_2_3, a_10_2_3); 00266 } 00267 catch (std::logic_error err) { 00268 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00269 *outStream << err.what() << '\n'; 00270 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00271 }; 00272 00273 *outStream << "-> absolute value with multidimensional arrays:\n"; 00274 00275 try { 00276 rst::absval(a_10_3, a_10_2_3); 00277 } 00278 catch (std::logic_error err) { 00279 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00280 *outStream << err.what() << '\n'; 00281 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00282 }; 00283 try { 00284 rst::absval(a_10_2_2, a_10_2_3); 00285 } 00286 catch (std::logic_error err) { 00287 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00288 *outStream << err.what() << '\n'; 00289 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00290 }; 00291 #endif 00292 00293 } 00294 catch (std::logic_error err) { 00295 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n"; 00296 *outStream << err.what() << '\n'; 00297 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00298 errorFlag = -1000; 00299 }; 00300 00301 00302 00303 *outStream \ 00304 << "\n" 00305 << "===============================================================================\n"\ 00306 << "| TEST 2: matrix / matrix-vector exceptions |\n"\ 00307 << "===============================================================================\n"; 00308 00309 try{ 00310 00311 FieldContainer<double> a_10_1_2_3_4(10, 1, 2, 3, 4); 00312 FieldContainer<double> b_10_1_2_3_4(10, 1, 2, 3, 4); 00313 FieldContainer<double> a_10(10); 00314 FieldContainer<double> a_9(9); 00315 FieldContainer<double> b_10(10); 00316 FieldContainer<double> a_10_15_4_4(10, 15, 4, 4); 00317 FieldContainer<double> b_10_15_4_4(10, 15, 4, 4); 00318 FieldContainer<double> a_10_2_2(10, 2, 2); 00319 FieldContainer<double> a_10_2_3(10, 2, 3); 00320 FieldContainer<double> b_10_2_3(10, 2, 3); 00321 00322 FieldContainer<double> a_1_1(1, 1); 00323 FieldContainer<double> b_1_1(1, 1); 00324 FieldContainer<double> a_2_2(2, 2); 00325 FieldContainer<double> b_2_2(2, 2); 00326 FieldContainer<double> a_3_3(3, 3); 00327 FieldContainer<double> b_3_3(3, 3); 00328 FieldContainer<double> a_2_3(2, 3); 00329 FieldContainer<double> a_4_4(4, 4); 00330 00331 FieldContainer<double> a_10_15_1_1(10, 15, 1, 1); 00332 FieldContainer<double> b_10_15_1_1(10, 15, 1, 1); 00333 FieldContainer<double> a_10_15_2_2(10, 15, 2, 2); 00334 FieldContainer<double> b_10_15_2_2(10, 15, 2, 2); 00335 FieldContainer<double> a_10_15_3_3(10, 15, 3, 3); 00336 FieldContainer<double> a_10_15_3_2(10, 15, 3, 2); 00337 FieldContainer<double> b_10_15_3_3(10, 15, 3, 3); 00338 FieldContainer<double> b_10_14(10, 14); 00339 FieldContainer<double> b_10_15(10, 15); 00340 FieldContainer<double> b_10_14_3(10, 14, 3); 00341 FieldContainer<double> b_10_15_3(10, 15, 3); 00342 00343 00344 #ifdef HAVE_INTREPID_DEBUG 00345 00346 *outStream << "-> inverse with multidimensional arrays:\n"; 00347 00348 try { 00349 rst::inverse(a_2_2, a_10_2_2); 00350 } 00351 catch (std::logic_error err) { 00352 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00353 *outStream << err.what() << '\n'; 00354 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00355 }; 00356 try { 00357 rst::inverse(b_10_1_2_3_4, a_10_1_2_3_4); 00358 } 00359 catch (std::logic_error err) { 00360 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00361 *outStream << err.what() << '\n'; 00362 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00363 }; 00364 try { 00365 rst::inverse(b_10, a_10); 00366 } 00367 catch (std::logic_error err) { 00368 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00369 *outStream << err.what() << '\n'; 00370 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00371 }; 00372 try { 00373 rst::inverse(a_10_2_2, a_10_2_3); 00374 } 00375 catch (std::logic_error err) { 00376 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00377 *outStream << err.what() << '\n'; 00378 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00379 }; 00380 try { 00381 rst::inverse(b_10_2_3, a_10_2_3); 00382 } 00383 catch (std::logic_error err) { 00384 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00385 *outStream << err.what() << '\n'; 00386 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00387 }; 00388 try { 00389 rst::inverse(b_10_15_4_4, a_10_15_4_4); 00390 } 00391 catch (std::logic_error err) { 00392 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00393 *outStream << err.what() << '\n'; 00394 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00395 }; 00396 try { 00397 rst::inverse(b_1_1, a_1_1); 00398 } 00399 catch (std::logic_error err) { 00400 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00401 *outStream << err.what() << '\n'; 00402 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00403 }; 00404 try { 00405 rst::inverse(b_2_2, a_2_2); 00406 } 00407 catch (std::logic_error err) { 00408 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00409 *outStream << err.what() << '\n'; 00410 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00411 }; 00412 try { 00413 rst::inverse(b_3_3, a_3_3); 00414 } 00415 catch (std::logic_error err) { 00416 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00417 *outStream << err.what() << '\n'; 00418 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00419 }; 00420 a_2_2[0] = 1.0; 00421 a_3_3[0] = 1.0; 00422 try { 00423 rst::inverse(b_2_2, a_2_2); 00424 } 00425 catch (std::logic_error err) { 00426 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00427 *outStream << err.what() << '\n'; 00428 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00429 }; 00430 try { 00431 rst::inverse(b_3_3, a_3_3); 00432 } 00433 catch (std::logic_error err) { 00434 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00435 *outStream << err.what() << '\n'; 00436 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00437 }; 00438 00439 *outStream << "-> transpose with multidimensional arrays:\n"; 00440 00441 try { 00442 rst::transpose(a_2_2, a_10_2_2); 00443 } 00444 catch (std::logic_error err) { 00445 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00446 *outStream << err.what() << '\n'; 00447 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00448 }; 00449 try { 00450 rst::transpose(b_10_1_2_3_4, a_10_1_2_3_4); 00451 } 00452 catch (std::logic_error err) { 00453 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00454 *outStream << err.what() << '\n'; 00455 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00456 }; 00457 try { 00458 rst::transpose(b_10, a_10); 00459 } 00460 catch (std::logic_error err) { 00461 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00462 *outStream << err.what() << '\n'; 00463 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00464 }; 00465 try { 00466 rst::transpose(a_10_2_2, a_10_2_3); 00467 } 00468 catch (std::logic_error err) { 00469 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00470 *outStream << err.what() << '\n'; 00471 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00472 }; 00473 try { 00474 rst::transpose(b_10_2_3, a_10_2_3); 00475 } 00476 catch (std::logic_error err) { 00477 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00478 *outStream << err.what() << '\n'; 00479 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00480 }; 00481 00482 *outStream << "-> determinant with multidimensional arrays:\n"; 00483 00484 try { 00485 rst::det(a_2_2, a_10_2_2); 00486 } 00487 catch (std::logic_error err) { 00488 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00489 *outStream << err.what() << '\n'; 00490 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00491 }; 00492 try { 00493 rst::det(a_10_2_2, a_10_1_2_3_4); 00494 } 00495 catch (std::logic_error err) { 00496 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00497 *outStream << err.what() << '\n'; 00498 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00499 }; 00500 try { 00501 rst::det(b_10_14, a_10_15_3_3); 00502 } 00503 catch (std::logic_error err) { 00504 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00505 *outStream << err.what() << '\n'; 00506 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00507 }; 00508 try { 00509 rst::det(a_9, a_10_2_2); 00510 } 00511 catch (std::logic_error err) { 00512 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00513 *outStream << err.what() << '\n'; 00514 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00515 }; 00516 try { 00517 rst::det(b_10, a_10_2_3); 00518 } 00519 catch (std::logic_error err) { 00520 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00521 *outStream << err.what() << '\n'; 00522 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00523 }; 00524 try { 00525 rst::det(b_10_15, a_10_15_4_4); 00526 } 00527 catch (std::logic_error err) { 00528 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00529 *outStream << err.what() << '\n'; 00530 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00531 }; 00532 try { 00533 rst::det(a_10_15_4_4); 00534 } 00535 catch (std::logic_error err) { 00536 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00537 *outStream << err.what() << '\n'; 00538 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00539 }; 00540 try { 00541 rst::det(a_2_3); 00542 } 00543 catch (std::logic_error err) { 00544 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00545 *outStream << err.what() << '\n'; 00546 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00547 }; 00548 try { 00549 rst::det(a_4_4); 00550 } 00551 catch (std::logic_error err) { 00552 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00553 *outStream << err.what() << '\n'; 00554 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00555 }; 00556 00557 *outStream << "-> matrix-vector product with multidimensional arrays:\n"; 00558 00559 try { 00560 rst::matvec(a_10_2_2, a_10_2_3, b_10_2_3); 00561 } 00562 catch (std::logic_error err) { 00563 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00564 *outStream << err.what() << '\n'; 00565 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00566 }; 00567 try { 00568 rst::matvec(a_2_2, a_2_2, a_10); 00569 } 00570 catch (std::logic_error err) { 00571 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00572 *outStream << err.what() << '\n'; 00573 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00574 }; 00575 try { 00576 rst::matvec(a_9, a_10_2_2, a_2_2); 00577 } 00578 catch (std::logic_error err) { 00579 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00580 *outStream << err.what() << '\n'; 00581 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00582 }; 00583 try { 00584 rst::matvec(b_10_15_3, a_10_15_3_3, b_10_14_3); 00585 } 00586 catch (std::logic_error err) { 00587 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00588 *outStream << err.what() << '\n'; 00589 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00590 }; 00591 try { 00592 rst::matvec(b_10_14_3, a_10_15_3_3, b_10_15_3); 00593 } 00594 catch (std::logic_error err) { 00595 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00596 *outStream << err.what() << '\n'; 00597 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00598 }; 00599 try { 00600 rst::matvec(b_10_15_3, a_10_15_3_2, b_10_15_3); 00601 } 00602 catch (std::logic_error err) { 00603 *outStream << "Expected Error ----------------------------------------------------------------\n"; 00604 *outStream << err.what() << '\n'; 00605 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00606 }; 00607 00608 #endif 00609 00610 } 00611 catch (std::logic_error err) { 00612 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n"; 00613 *outStream << err.what() << '\n'; 00614 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00615 errorFlag = -1000; 00616 }; 00617 00618 #ifdef HAVE_INTREPID_DEBUG 00619 if (Teuchos::TestForException_getThrowNumber() != endThrowNumber) 00620 errorFlag++; 00621 #endif 00622 00623 00624 *outStream \ 00625 << "\n" 00626 << "===============================================================================\n"\ 00627 << "| TEST 2: correctness of math operations |\n"\ 00628 << "===============================================================================\n"; 00629 00630 outStream->precision(20); 00631 00632 try{ 00633 00634 // two-dimensional base containers 00635 for (int dim=3; dim>0; dim--) { 00636 int i0=4, i1=5; 00637 FieldContainer<double> ma_x_x_d_d(i0, i1, dim, dim); 00638 FieldContainer<double> mb_x_x_d_d(i0, i1, dim, dim); 00639 FieldContainer<double> mc_x_x_d_d(i0, i1, dim, dim); 00640 FieldContainer<double> va_x_x_d(i0, i1, dim); 00641 FieldContainer<double> vb_x_x_d(i0, i1, dim); 00642 FieldContainer<double> vc_x_x_d(i0, i1, dim); 00643 FieldContainer<double> vdot_x_x(i0, i1); 00644 FieldContainer<double> vnorms_x_x(i0, i1); 00645 FieldContainer<double> vnorms_x(i0); 00646 double zero = INTREPID_TOL*100.0; 00647 00648 // fill with random numbers 00649 for (int i=0; i<ma_x_x_d_d.size(); i++) { 00650 ma_x_x_d_d[i] = Teuchos::ScalarTraits<double>::random(); 00651 } 00652 for (int i=0; i<va_x_x_d.size(); i++) { 00653 va_x_x_d[i] = Teuchos::ScalarTraits<double>::random(); 00654 } 00655 00656 *outStream << "\n************ Checking vectorNorm ************\n"; 00657 00658 rst::vectorNorm(vnorms_x_x, va_x_x_d, NORM_TWO); 00659 *outStream << va_x_x_d; 00660 *outStream << vnorms_x_x; 00661 if ( std::abs(rst::vectorNorm(&vnorms_x_x[0], vnorms_x_x.size(), NORM_TWO) - 00662 rst::vectorNorm(&va_x_x_d[0], va_x_x_d.size(), NORM_TWO)) > zero) { 00663 *outStream << "\n\nINCORRECT vectorNorm NORM_TWO\n\n"; 00664 errorFlag = -1000; 00665 } 00666 00667 rst::vectorNorm(vnorms_x_x, va_x_x_d, NORM_ONE); 00668 *outStream << va_x_x_d; 00669 *outStream << vnorms_x_x; 00670 if ( std::abs(rst::vectorNorm(&vnorms_x_x[0], vnorms_x_x.size(), NORM_ONE) - 00671 rst::vectorNorm(&va_x_x_d[0], va_x_x_d.size(), NORM_ONE)) > zero) { 00672 *outStream << "\n\nINCORRECT vectorNorm NORM_ONE\n\n"; 00673 errorFlag = -1000; 00674 } 00675 00676 rst::vectorNorm(vnorms_x_x, va_x_x_d, NORM_INF); 00677 *outStream << va_x_x_d; 00678 *outStream << vnorms_x_x; 00679 if ( std::abs(rst::vectorNorm(&vnorms_x_x[0], vnorms_x_x.size(), NORM_INF) - 00680 rst::vectorNorm(&va_x_x_d[0], va_x_x_d.size(), NORM_INF)) > zero) { 00681 *outStream << "\n\nINCORRECT vectorNorm NORM_INF\n\n"; 00682 errorFlag = -1000; 00683 } 00684 00685 /******************************************/ 00686 00687 00688 *outStream << "\n************ Checking inverse, subtract, and vectorNorm ************\n"; 00689 00690 rst::inverse(mb_x_x_d_d, ma_x_x_d_d); // B = inv(A) 00691 rst::inverse(mc_x_x_d_d, mb_x_x_d_d); // C = inv(B) ~= A 00692 *outStream << ma_x_x_d_d << mb_x_x_d_d << mc_x_x_d_d; 00693 00694 rst::subtract(&mc_x_x_d_d[0], &ma_x_x_d_d[0], ma_x_x_d_d.size()); // C = C - A ~= 0 00695 00696 if (rst::vectorNorm(&mc_x_x_d_d[0], mc_x_x_d_d.size(), NORM_ONE) > zero) { 00697 *outStream << "\n\nINCORRECT inverse OR subtract OR vectorNorm\n\n"; 00698 errorFlag = -1000; 00699 } 00700 00701 /******************************************/ 00702 00703 00704 *outStream << "\n********** Checking determinant **********\n"; 00705 00706 FieldContainer<double> detA_x_x(i0, i1); 00707 FieldContainer<double> detB_x_x(i0, i1); 00708 00709 rst::det(detA_x_x, ma_x_x_d_d); 00710 rst::det(detB_x_x, mb_x_x_d_d); 00711 *outStream << detA_x_x << detB_x_x; 00712 00713 if ( (rst::dot(&detA_x_x[0], &detB_x_x[0], detA_x_x.size()) - (double)(i0*i1)) > zero) { 00714 *outStream << "\n\nINCORRECT det\n\n" ; 00715 errorFlag = -1000; 00716 } 00717 00718 *outStream << "\n det(A)*det(inv(A)) = " << 00719 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)) 00720 << "\n"; 00721 00722 if ( (rst::det(&ma_x_x_d_d[0], ma_x_x_d_d.dimension(3))* 00723 rst::det(&mb_x_x_d_d[0], mb_x_x_d_d.dimension(3)) - (double)1) > zero) { 00724 *outStream << "\n\nINCORRECT det\n\n" ; 00725 errorFlag = -1000; 00726 } 00727 00728 /******************************************/ 00729 00730 00731 *outStream << "\n************ Checking transpose and subtract ************\n"; 00732 00733 rst::transpose(mb_x_x_d_d, ma_x_x_d_d); // B = A^T 00734 rst::transpose(mc_x_x_d_d, mb_x_x_d_d); // C = B^T = A 00735 *outStream << ma_x_x_d_d << mb_x_x_d_d << mc_x_x_d_d; 00736 00737 rst::subtract(&mc_x_x_d_d[0], &ma_x_x_d_d[0], ma_x_x_d_d.size()); // C = C - A = 0 00738 00739 if (rst::vectorNorm(&mc_x_x_d_d[0], mc_x_x_d_d.size(), NORM_ONE) > zero) { 00740 *outStream << "\n\nINCORRECT transpose OR subtract OR vectorNorm\n\n" ; 00741 errorFlag = -1000; 00742 } 00743 00744 /******************************************/ 00745 00746 00747 *outStream << "\n************ Checking matvec, vectorNorm, subtract, and inverse ************\n"; 00748 00749 rst::inverse(mb_x_x_d_d, ma_x_x_d_d); // B = inv(A) 00750 rst::inverse(mc_x_x_d_d, mb_x_x_d_d); // C = inv(B) ~= A 00751 rst::matvec(vb_x_x_d, ma_x_x_d_d, va_x_x_d); // b = A*a 00752 rst::matvec(vc_x_x_d, mb_x_x_d_d, vb_x_x_d); // c = inv(A)*(A*a) ~= a 00753 rst::subtract(vc_x_x_d, va_x_x_d); // c = c - a ~= 0 00754 *outStream << vc_x_x_d; 00755 00756 rst::vectorNorm(vnorms_x_x, vc_x_x_d, NORM_ONE); 00757 rst::vectorNorm(vnorms_x, vnorms_x_x, NORM_INF); 00758 if (rst::vectorNorm(vnorms_x, NORM_TWO) > zero) { 00759 *outStream << "\n\nINCORRECT matvec OR inverse OR subtract OR vectorNorm\n\n"; 00760 errorFlag = -1000; 00761 } 00762 00763 /******************************************/ 00764 00765 00766 *outStream << "\n************ Checking add, subtract, absval, and scale ************\n"; 00767 00768 double x = 1.234; 00769 rst::add(vc_x_x_d, va_x_x_d, vb_x_x_d); // c = a + b 00770 rst::subtract(vc_x_x_d, vb_x_x_d); // c = c - b = a 00771 rst::scale(vb_x_x_d, vc_x_x_d, x); // b = c*x; 00772 rst::scale(vc_x_x_d, vb_x_x_d, (1.0/x)); // c = b*(1/x) = a; 00773 rst::subtract(vb_x_x_d, vc_x_x_d, va_x_x_d); // b = c - a ~= 0 00774 rst::absval(vc_x_x_d, vb_x_x_d); // c = |b| 00775 rst::scale(vb_x_x_d, vc_x_x_d, -1.0); // b = -c 00776 rst::absval(vc_x_x_d, vb_x_x_d); // c = |b| 00777 rst::add(vc_x_x_d, vb_x_x_d); // c = c + b === 0 00778 *outStream << vc_x_x_d; 00779 00780 rst::vectorNorm(vnorms_x_x, vc_x_x_d, NORM_ONE); 00781 rst::vectorNorm(vnorms_x, vnorms_x_x, NORM_INF); 00782 if (rst::vectorNorm(vnorms_x, NORM_TWO) > (double)0) { 00783 *outStream << "\n\nSign flips combined with std::abs might not be invertible on this platform!\n" 00784 << "Potential IEEE compliance issues!\n\n"; 00785 if (rst::vectorNorm(vnorms_x, NORM_TWO) > zero) { 00786 *outStream << "\n\nINCORRECT add OR subtract OR scale OR absval OR vectorNorm\n\n"; 00787 errorFlag = -1000; 00788 } 00789 } 00790 00791 /******************************************/ 00792 00793 00794 *outStream << "\n************ Checking dot and vectorNorm ************\n"; 00795 00796 for (int i=0; i<va_x_x_d.size(); i++) { 00797 va_x_x_d[i] = 2.0; 00798 } 00799 rst::dot(vdot_x_x, va_x_x_d, va_x_x_d); // dot = a'*a 00800 *outStream << vdot_x_x; 00801 00802 rst::vectorNorm(vnorms_x, vdot_x_x, NORM_ONE); 00803 if (rst::vectorNorm(vnorms_x, NORM_ONE) - (double)(4.0*dim*i0*i1) > zero) { 00804 *outStream << "\n\nINCORRECT dot OR vectorNorm\n\n"; 00805 errorFlag = -1000; 00806 } 00807 00808 /******************************************/ 00809 00810 *outStream << "\n"; 00811 } 00812 00813 // one-dimensional base containers 00814 for (int dim=3; dim>0; dim--) { 00815 int i0=7; 00816 FieldContainer<double> ma_x_d_d(i0, dim, dim); 00817 FieldContainer<double> mb_x_d_d(i0, dim, dim); 00818 FieldContainer<double> mc_x_d_d(i0, dim, dim); 00819 FieldContainer<double> va_x_d(i0, dim); 00820 FieldContainer<double> vb_x_d(i0, dim); 00821 FieldContainer<double> vc_x_d(i0, dim); 00822 FieldContainer<double> vdot_x(i0); 00823 FieldContainer<double> vnorms_x(i0); 00824 double zero = INTREPID_TOL*100.0; 00825 00826 // fill with random numbers 00827 for (int i=0; i<ma_x_d_d.size(); i++) { 00828 ma_x_d_d[i] = Teuchos::ScalarTraits<double>::random(); 00829 } 00830 for (int i=0; i<va_x_d.size(); i++) { 00831 va_x_d[i] = Teuchos::ScalarTraits<double>::random(); 00832 } 00833 00834 *outStream << "\n************ Checking vectorNorm ************\n"; 00835 00836 rst::vectorNorm(vnorms_x, va_x_d, NORM_TWO); 00837 *outStream << va_x_d; 00838 *outStream << vnorms_x; 00839 if ( std::abs(rst::vectorNorm(&vnorms_x[0], vnorms_x.size(), NORM_TWO) - 00840 rst::vectorNorm(&va_x_d[0], va_x_d.size(), NORM_TWO)) > zero) { 00841 *outStream << "\n\nINCORRECT vectorNorm NORM_TWO\n\n"; 00842 errorFlag = -1000; 00843 } 00844 00845 rst::vectorNorm(vnorms_x, va_x_d, NORM_ONE); 00846 *outStream << va_x_d; 00847 *outStream << vnorms_x; 00848 if ( std::abs(rst::vectorNorm(&vnorms_x[0], vnorms_x.size(), NORM_ONE) - 00849 rst::vectorNorm(&va_x_d[0], va_x_d.size(), NORM_ONE)) > zero) { 00850 *outStream << "\n\nINCORRECT vectorNorm NORM_ONE\n\n"; 00851 errorFlag = -1000; 00852 } 00853 00854 rst::vectorNorm(vnorms_x, va_x_d, NORM_INF); 00855 *outStream << va_x_d; 00856 *outStream << vnorms_x; 00857 if ( std::abs(rst::vectorNorm(&vnorms_x[0], vnorms_x.size(), NORM_INF) - 00858 rst::vectorNorm(&va_x_d[0], va_x_d.size(), NORM_INF)) > zero) { 00859 *outStream << "\n\nINCORRECT vectorNorm NORM_INF\n\n"; 00860 errorFlag = -1000; 00861 } 00862 00863 /******************************************/ 00864 00865 00866 *outStream << "\n************ Checking inverse, subtract, and vectorNorm ************\n"; 00867 00868 rst::inverse(mb_x_d_d, ma_x_d_d); // B = inv(A) 00869 rst::inverse(mc_x_d_d, mb_x_d_d); // C = inv(B) ~= A 00870 *outStream << ma_x_d_d << mb_x_d_d << mc_x_d_d; 00871 00872 rst::subtract(&mc_x_d_d[0], &ma_x_d_d[0], ma_x_d_d.size()); // C = C - A ~= 0 00873 00874 if (rst::vectorNorm(&mc_x_d_d[0], mc_x_d_d.size(), NORM_ONE) > zero) { 00875 *outStream << "\n\nINCORRECT inverse OR subtract OR vectorNorm\n\n"; 00876 errorFlag = -1000; 00877 } 00878 00879 /******************************************/ 00880 00881 00882 *outStream << "\n********** Checking determinant **********\n"; 00883 00884 FieldContainer<double> detA_x(i0); 00885 FieldContainer<double> detB_x(i0); 00886 00887 rst::det(detA_x, ma_x_d_d); 00888 rst::det(detB_x, mb_x_d_d); 00889 *outStream << detA_x << detB_x; 00890 00891 if ( (rst::dot(&detA_x[0], &detB_x[0], detA_x.size()) - (double)i0) > zero) { 00892 *outStream << "\n\nINCORRECT det\n\n" ; 00893 errorFlag = -1000; 00894 } 00895 00896 *outStream << "\n det(A)*det(inv(A)) = " << 00897 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)) 00898 << "\n"; 00899 00900 if ( (rst::det(&ma_x_d_d[0], ma_x_d_d.dimension(2))* 00901 rst::det(&mb_x_d_d[0], mb_x_d_d.dimension(2)) - (double)1) > zero) { 00902 *outStream << "\n\nINCORRECT det\n\n" ; 00903 errorFlag = -1000; 00904 } 00905 00906 /******************************************/ 00907 00908 00909 *outStream << "\n************ Checking transpose and subtract ************\n"; 00910 00911 rst::transpose(mb_x_d_d, ma_x_d_d); // B = A^T 00912 rst::transpose(mc_x_d_d, mb_x_d_d); // C = B^T = A 00913 *outStream << ma_x_d_d << mb_x_d_d << mc_x_d_d; 00914 00915 rst::subtract(&mc_x_d_d[0], &ma_x_d_d[0], ma_x_d_d.size()); // C = C - A = 0 00916 00917 if (rst::vectorNorm(&mc_x_d_d[0], mc_x_d_d.size(), NORM_ONE) > zero) { 00918 *outStream << "\n\nINCORRECT transpose OR subtract OR vectorNorm\n\n" ; 00919 errorFlag = -1000; 00920 } 00921 00922 /******************************************/ 00923 00924 00925 *outStream << "\n************ Checking matvec, vectorNorm, subtract, and inverse ************\n"; 00926 00927 rst::inverse(mb_x_d_d, ma_x_d_d); // B = inv(A) 00928 rst::inverse(mc_x_d_d, mb_x_d_d); // C = inv(B) ~= A 00929 rst::matvec(vb_x_d, ma_x_d_d, va_x_d); // b = A*a 00930 rst::matvec(vc_x_d, mb_x_d_d, vb_x_d); // c = inv(A)*(A*a) ~= a 00931 rst::subtract(vc_x_d, va_x_d); // c = c - a ~= 0 00932 *outStream << vc_x_d; 00933 00934 rst::vectorNorm(vnorms_x, vc_x_d, NORM_ONE); 00935 if (rst::vectorNorm(vnorms_x, NORM_TWO) > zero) { 00936 *outStream << "\n\nINCORRECT matvec OR inverse OR subtract OR vectorNorm\n\n"; 00937 errorFlag = -1000; 00938 } 00939 00940 /******************************************/ 00941 00942 00943 *outStream << "\n************ Checking add, subtract, absval, and scale ************\n"; 00944 00945 double x = 1.234; 00946 rst::add(vc_x_d, va_x_d, vb_x_d); // c = a + b 00947 rst::subtract(vc_x_d, vb_x_d); // c = c - b = a 00948 rst::scale(vb_x_d, vc_x_d, x); // b = c*x; 00949 rst::scale(vc_x_d, vb_x_d, (1.0/x)); // c = b*(1/x) = a; 00950 rst::subtract(vb_x_d, vc_x_d, va_x_d); // b = c - a ~= 0 00951 rst::absval(vc_x_d, vb_x_d); // c = |b| 00952 rst::scale(vb_x_d, vc_x_d, -1.0); // b = -c 00953 rst::absval(vc_x_d, vb_x_d); // c = |b| 00954 rst::add(vc_x_d, vb_x_d); // c = c + b === 0 00955 *outStream << vc_x_d; 00956 00957 rst::vectorNorm(vnorms_x, vc_x_d, NORM_ONE); 00958 if (rst::vectorNorm(vnorms_x, NORM_TWO) > (double)0) { 00959 *outStream << "\n\nSign flips combined with std::abs might not be invertible on this platform!\n" 00960 << "Potential IEEE compliance issues!\n\n"; 00961 if (rst::vectorNorm(vnorms_x, NORM_TWO) > zero) { 00962 *outStream << "\n\nINCORRECT add OR subtract OR scale OR absval OR vectorNorm\n\n"; 00963 errorFlag = -1000; 00964 } 00965 } 00966 00967 /******************************************/ 00968 00969 00970 *outStream << "\n************ Checking dot and vectorNorm ************\n"; 00971 00972 for (int i=0; i<va_x_d.size(); i++) { 00973 va_x_d[i] = 2.0; 00974 } 00975 rst::dot(vdot_x, va_x_d, va_x_d); // dot = a'*a 00976 *outStream << vdot_x; 00977 00978 if (rst::vectorNorm(vdot_x, NORM_ONE) - (double)(4.0*dim*i0) > zero) { 00979 *outStream << "\n\nINCORRECT dot OR vectorNorm\n\n"; 00980 errorFlag = -1000; 00981 } 00982 00983 /******************************************/ 00984 00985 *outStream << "\n"; 00986 } 00987 } 00988 catch (std::logic_error err) { 00989 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n"; 00990 *outStream << err.what() << '\n'; 00991 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00992 errorFlag = -1000; 00993 }; 00994 00995 if (errorFlag != 0) 00996 std::cout << "End Result: TEST FAILED\n"; 00997 else 00998 std::cout << "End Result: TEST PASSED\n"; 00999 01000 // reset format state of std::cout 01001 std::cout.copyfmt(oldFormatState); 01002 01003 return errorFlag; 01004 }
1.7.6.1