Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/test/Shared/RealSpaceTools/test_01.cpp
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 }