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