Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/test/Shared/ArrayTools/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_ArrayTools.hpp"
00050 #include "Intrepid_FieldContainer.hpp"
00051 #include "Intrepid_RealSpaceTools.hpp"
00052 #include "Teuchos_oblackholestream.hpp"
00053 #include "Teuchos_RCP.hpp"
00054 #include "Teuchos_ScalarTraits.hpp"
00055 #include "Teuchos_GlobalMPISession.hpp"
00056 
00057 using namespace std;
00058 using namespace Intrepid;
00059 
00060 #define INTREPID_TEST_COMMAND( S )                                                                                  \
00061 {                                                                                                                   \
00062   try {                                                                                                             \
00063     S ;                                                                                                             \
00064   }                                                                                                                 \
00065   catch (std::logic_error err) {                                                                                    \
00066       *outStream << "Expected Error ----------------------------------------------------------------\n";            \
00067       *outStream << err.what() << '\n';                                                                             \
00068       *outStream << "-------------------------------------------------------------------------------" << "\n\n";    \
00069   };                                                                                                                \
00070 }
00071 
00072 
00073 int main(int argc, char *argv[]) {
00074 
00075   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00076 
00077   // This little trick lets us print to std::cout only if
00078   // a (dummy) command-line argument is provided.
00079   int iprint     = argc - 1;
00080   Teuchos::RCP<std::ostream> outStream;
00081   Teuchos::oblackholestream bhs; // outputs nothing
00082   if (iprint > 0)
00083     outStream = Teuchos::rcp(&std::cout, false);
00084   else
00085     outStream = Teuchos::rcp(&bhs, false);
00086 
00087   // Save the format state of the original std::cout.
00088   Teuchos::oblackholestream oldFormatState;
00089   oldFormatState.copyfmt(std::cout);
00090 
00091   *outStream \
00092   << "===============================================================================\n" \
00093   << "|                                                                             |\n" \
00094   << "|                       Unit Test (ArrayTools)                                |\n" \
00095   << "|                                                                             |\n" \
00096   << "|     1) Array operations: contractions                                       |\n" \
00097   << "|                                                                             |\n" \
00098   << "|  Questions? Contact  Pavel Bochev (pbboche@sandia.gov) or                   |\n" \
00099   << "|                      Denis Ridzal (dridzal@sandia.gov).                     |\n" \
00100   << "|                                                                             |\n" \
00101   << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00102   << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00103   << "|                                                                             |\n" \
00104   << "===============================================================================\n";
00105 
00106 
00107   int errorFlag = 0;
00108 #ifdef HAVE_INTREPID_DEBUG
00109   int beginThrowNumber = Teuchos::TestForException_getThrowNumber();
00110   int endThrowNumber = beginThrowNumber + 81;
00111 #endif
00112 
00113   typedef ArrayTools art; 
00114   typedef RealSpaceTools<double> rst; 
00115 #ifdef HAVE_INTREPID_DEBUG
00116   ArrayTools atools;
00117 #endif
00118 
00119   *outStream \
00120   << "\n"
00121   << "===============================================================================\n"\
00122   << "| TEST 1: exceptions                                                          |\n"\
00123   << "===============================================================================\n";
00124 
00125   try{
00126 
00127 #ifdef HAVE_INTREPID_DEBUG
00128     FieldContainer<double> a_2_2(2, 2);
00129     FieldContainer<double> a_10_2(10, 2);
00130     FieldContainer<double> a_10_3(10, 3);
00131     FieldContainer<double> a_10_2_2(10, 2, 2);
00132     FieldContainer<double> a_10_2_3(10, 2, 3);
00133     FieldContainer<double> a_10_3_2(10, 3, 2);
00134     FieldContainer<double> a_9_2_2(9, 2, 2);
00135 
00136     *outStream << "-> contractFieldFieldScalar:\n";
00137     INTREPID_TEST_COMMAND( atools.contractFieldFieldScalar<double>(a_2_2, a_2_2, a_2_2, COMP_CPP) );
00138     INTREPID_TEST_COMMAND( atools.contractFieldFieldScalar<double>(a_2_2, a_10_2_2, a_2_2, COMP_CPP) );
00139     INTREPID_TEST_COMMAND( atools.contractFieldFieldScalar<double>(a_2_2, a_10_2_2, a_10_2_2, COMP_CPP) );
00140     INTREPID_TEST_COMMAND( atools.contractFieldFieldScalar<double>(a_10_2_2, a_9_2_2, a_10_2_2, COMP_CPP) );
00141     INTREPID_TEST_COMMAND( atools.contractFieldFieldScalar<double>(a_10_2_2, a_10_2_2, a_10_2_3, COMP_CPP) );
00142     INTREPID_TEST_COMMAND( atools.contractFieldFieldScalar<double>(a_9_2_2, a_10_3_2, a_10_2_2, COMP_CPP) );
00143     INTREPID_TEST_COMMAND( atools.contractFieldFieldScalar<double>(a_10_3_2, a_10_2_2, a_10_3_2, COMP_CPP) );
00144     INTREPID_TEST_COMMAND( atools.contractFieldFieldScalar<double>(a_10_2_2, a_10_2_2, a_10_3_2, COMP_CPP) );
00145     INTREPID_TEST_COMMAND( atools.contractFieldFieldScalar<double>(a_10_2_3, a_10_2_2, a_10_3_2, COMP_ENGINE_MAX) );
00146     INTREPID_TEST_COMMAND( atools.contractFieldFieldScalar<double>(a_10_2_3, a_10_2_2, a_10_3_2, COMP_CPP) );
00147 
00148     FieldContainer<double> a_10_2_2_2(10, 2, 2, 2);
00149     FieldContainer<double> a_9_2_2_2(9, 2, 2, 2);
00150     FieldContainer<double> a_10_3_2_2(10, 3, 2, 2);
00151     FieldContainer<double> a_10_2_3_2(10, 2, 3, 2);
00152     FieldContainer<double> a_10_2_2_3(10, 2, 2, 3);
00153 
00154     *outStream << "-> contractFieldFieldVector:\n";
00155     INTREPID_TEST_COMMAND( atools.contractFieldFieldVector<double>(a_2_2, a_2_2, a_2_2, COMP_CPP) );
00156     INTREPID_TEST_COMMAND( atools.contractFieldFieldVector<double>(a_2_2, a_10_2_2_2, a_2_2, COMP_CPP) );
00157     INTREPID_TEST_COMMAND( atools.contractFieldFieldVector<double>(a_2_2, a_10_2_2_2, a_10_2_2_2, COMP_CPP) );
00158     INTREPID_TEST_COMMAND( atools.contractFieldFieldVector<double>(a_10_2_2, a_9_2_2_2, a_10_2_2_2, COMP_CPP) );
00159     INTREPID_TEST_COMMAND( atools.contractFieldFieldVector<double>(a_10_2_2, a_10_2_2_2, a_10_2_3_2, COMP_CPP) );
00160     INTREPID_TEST_COMMAND( atools.contractFieldFieldVector<double>(a_10_2_2, a_10_2_2_2, a_10_2_2_3, COMP_CPP) );
00161     INTREPID_TEST_COMMAND( atools.contractFieldFieldVector<double>(a_9_2_2, a_10_2_2_2, a_10_2_2_2, COMP_CPP) );
00162     INTREPID_TEST_COMMAND( atools.contractFieldFieldVector<double>(a_10_3_2, a_10_2_2_2, a_10_3_2_2, COMP_CPP) );
00163     INTREPID_TEST_COMMAND( atools.contractFieldFieldVector<double>(a_10_2_2, a_10_2_2_2, a_10_3_2_2, COMP_CPP) );
00164     INTREPID_TEST_COMMAND( atools.contractFieldFieldVector<double>(a_10_2_3, a_10_2_2_2, a_10_3_2_2, COMP_ENGINE_MAX) );
00165     INTREPID_TEST_COMMAND( atools.contractFieldFieldVector<double>(a_10_2_3, a_10_2_2_2, a_10_3_2_2, COMP_CPP) );
00166 
00167     FieldContainer<double> a_10_2_2_2_2(10, 2, 2, 2, 2);
00168     FieldContainer<double> a_9_2_2_2_2(9, 2, 2, 2, 2);
00169     FieldContainer<double> a_10_3_2_2_2(10, 3, 2, 2, 2);
00170     FieldContainer<double> a_10_2_3_2_2(10, 2, 3, 2, 2);
00171     FieldContainer<double> a_10_2_2_3_2(10, 2, 2, 3, 2);
00172     FieldContainer<double> a_10_2_2_2_3(10, 2, 2, 2, 3);
00173 
00174     *outStream << "-> contractFieldFieldTensor:\n";
00175     INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<double>(a_2_2, a_2_2, a_2_2, COMP_CPP) );
00176     INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<double>(a_2_2, a_10_2_2_2_2, a_2_2, COMP_CPP) );
00177     INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<double>(a_2_2, a_10_2_2_2_2, a_10_2_2_2_2, COMP_CPP) );
00178     INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<double>(a_10_2_2, a_9_2_2_2_2, a_10_2_2_2_2, COMP_CPP) );
00179     INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<double>(a_10_2_2, a_10_2_2_2_2, a_10_2_3_2_2, COMP_CPP) );
00180     INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<double>(a_10_2_2, a_10_2_2_2_2, a_10_2_2_3_2, COMP_CPP) );
00181     INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<double>(a_10_2_2, a_10_2_2_2_2, a_10_2_2_2_3, COMP_CPP) );
00182     INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<double>(a_9_2_2, a_10_2_2_2_2, a_10_2_2_2_2, COMP_CPP) );
00183     INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<double>(a_10_3_2, a_10_2_2_2_2, a_10_2_2_2_2, COMP_CPP) );
00184     INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<double>(a_10_2_2, a_10_2_2_2_2, a_10_3_2_2_2, COMP_CPP) );
00185     INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<double>(a_10_2_3, a_10_2_2_2_2, a_10_3_2_2_2, COMP_ENGINE_MAX) );
00186     INTREPID_TEST_COMMAND( atools.contractFieldFieldTensor<double>(a_10_2_3, a_10_2_2_2_2, a_10_3_2_2_2, COMP_CPP) );
00187 
00188     FieldContainer<double> a_9_2(9, 2);
00189     FieldContainer<double> a_10_1(10, 1);
00190 
00191     *outStream << "-> contractDataFieldScalar:\n";
00192     INTREPID_TEST_COMMAND( atools.contractDataFieldScalar<double>(a_2_2, a_2_2, a_2_2, COMP_CPP) );
00193     INTREPID_TEST_COMMAND( atools.contractDataFieldScalar<double>(a_2_2, a_10_2_2, a_10_2_2, COMP_CPP) );
00194     INTREPID_TEST_COMMAND( atools.contractDataFieldScalar<double>(a_10_2_2, a_2_2, a_10_2_2, COMP_CPP) );
00195     INTREPID_TEST_COMMAND( atools.contractDataFieldScalar<double>(a_2_2, a_10_2, a_9_2_2, COMP_CPP) );
00196     INTREPID_TEST_COMMAND( atools.contractDataFieldScalar<double>(a_2_2, a_10_3, a_10_2_2, COMP_CPP) );
00197     INTREPID_TEST_COMMAND( atools.contractDataFieldScalar<double>(a_9_2, a_10_2, a_10_2_2, COMP_CPP) );
00198     INTREPID_TEST_COMMAND( atools.contractDataFieldScalar<double>(a_10_2, a_10_2, a_10_3_2, COMP_CPP) );
00199     INTREPID_TEST_COMMAND( atools.contractDataFieldScalar<double>(a_10_2, a_10_2, a_10_2_2, COMP_ENGINE_MAX) );
00200     INTREPID_TEST_COMMAND( atools.contractDataFieldScalar<double>(a_10_2, a_10_2, a_10_2_2, COMP_CPP) );
00201     INTREPID_TEST_COMMAND( atools.contractDataFieldScalar<double>(a_10_2, a_10_1, a_10_2_2, COMP_CPP) );
00202 
00203     FieldContainer<double> a_10_1_2(10, 1, 2);
00204     FieldContainer<double> a_10_1_3(10, 1, 3);
00205 
00206     *outStream << "-> contractDataFieldVector:\n";
00207     INTREPID_TEST_COMMAND( atools.contractDataFieldVector<double>(a_2_2, a_2_2, a_2_2, COMP_CPP) );
00208     INTREPID_TEST_COMMAND( atools.contractDataFieldVector<double>(a_2_2, a_2_2, a_10_2_2_2, COMP_CPP) );
00209     INTREPID_TEST_COMMAND( atools.contractDataFieldVector<double>(a_10_2_2, a_10_2_2, a_10_2_2_2, COMP_CPP) );
00210     INTREPID_TEST_COMMAND( atools.contractDataFieldVector<double>(a_2_2, a_10_2_2, a_9_2_2_2, COMP_CPP) );
00211     INTREPID_TEST_COMMAND( atools.contractDataFieldVector<double>(a_2_2, a_10_2_2, a_10_2_3_2, COMP_CPP) );
00212     INTREPID_TEST_COMMAND( atools.contractDataFieldVector<double>(a_2_2, a_10_2_2, a_10_2_2_3, COMP_CPP) );
00213     INTREPID_TEST_COMMAND( atools.contractDataFieldVector<double>(a_9_2, a_10_2_2, a_10_2_2_2, COMP_CPP) );
00214     INTREPID_TEST_COMMAND( atools.contractDataFieldVector<double>(a_10_2, a_10_2_2, a_10_3_2_2, COMP_CPP) );
00215     INTREPID_TEST_COMMAND( atools.contractDataFieldVector<double>(a_10_2, a_10_2_2, a_10_2_2_2, COMP_ENGINE_MAX) );
00216     INTREPID_TEST_COMMAND( atools.contractDataFieldVector<double>(a_10_2, a_10_2_2, a_10_2_2_2, COMP_CPP) );
00217     INTREPID_TEST_COMMAND( atools.contractDataFieldVector<double>(a_10_2, a_10_1_2, a_10_2_2_2, COMP_CPP) );
00218 
00219     FieldContainer<double> a_10_1_2_2(10, 1, 2, 2);
00220 
00221     *outStream << "-> contractDataFieldTensor:\n";
00222     INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<double>(a_2_2, a_2_2, a_2_2, COMP_CPP) );
00223     INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<double>(a_2_2, a_2_2, a_10_2_2_2_2, COMP_CPP) );
00224     INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<double>(a_10_2_2, a_10_2_2_2, a_10_2_2_2_2, COMP_CPP) );
00225     INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<double>(a_2_2, a_10_2_2_2, a_9_2_2_2_2, COMP_CPP) );
00226     INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<double>(a_2_2, a_10_2_2_2, a_10_2_3_2_2, COMP_CPP) );
00227     INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<double>(a_2_2, a_10_2_2_2, a_10_2_2_3_2, COMP_CPP) );
00228     INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<double>(a_2_2, a_10_2_2_2, a_10_2_2_2_3, COMP_CPP) );
00229     INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<double>(a_9_2, a_10_2_2_2, a_10_2_2_2_2, COMP_CPP) );
00230     INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<double>(a_10_2, a_10_2_2_2, a_10_3_2_2_2, COMP_CPP) );
00231     INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<double>(a_10_2, a_10_2_2_2, a_10_2_2_2_2, COMP_ENGINE_MAX) );
00232     INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<double>(a_10_2, a_10_2_2_2, a_10_2_2_2_2, COMP_CPP) );
00233     INTREPID_TEST_COMMAND( atools.contractDataFieldTensor<double>(a_10_2, a_10_1_2_2, a_10_2_2_2_2, COMP_CPP) );
00234 
00235     FieldContainer<double> a_2(2);
00236     FieldContainer<double> a_10(10);
00237 
00238     *outStream << "-> contractDataDataScalar:\n";
00239     INTREPID_TEST_COMMAND( atools.contractDataDataScalar<double>(a_2_2, a_10_2_2, a_10_2_2, COMP_CPP) );
00240     INTREPID_TEST_COMMAND( atools.contractDataDataScalar<double>(a_2_2, a_10_2, a_10_2_2, COMP_CPP) );
00241     INTREPID_TEST_COMMAND( atools.contractDataDataScalar<double>(a_2_2, a_10_2, a_10_2, COMP_CPP) );
00242     INTREPID_TEST_COMMAND( atools.contractDataDataScalar<double>(a_2, a_9_2, a_10_2, COMP_CPP) );
00243     INTREPID_TEST_COMMAND( atools.contractDataDataScalar<double>(a_2, a_10_2, a_10_3, COMP_CPP) );
00244     INTREPID_TEST_COMMAND( atools.contractDataDataScalar<double>(a_2, a_10_2, a_10_2, COMP_CPP) );
00245     INTREPID_TEST_COMMAND( atools.contractDataDataScalar<double>(a_10, a_10_2, a_10_2, COMP_ENGINE_MAX) );
00246     INTREPID_TEST_COMMAND( atools.contractDataDataScalar<double>(a_10, a_10_2, a_10_2, COMP_CPP) );
00247 
00248     *outStream << "-> contractDataDataVector:\n";
00249     INTREPID_TEST_COMMAND( atools.contractDataDataVector<double>(a_2_2, a_2_2, a_2_2, COMP_CPP) );
00250     INTREPID_TEST_COMMAND( atools.contractDataDataVector<double>(a_2_2, a_10_2_2, a_2_2, COMP_CPP) );
00251     INTREPID_TEST_COMMAND( atools.contractDataDataVector<double>(a_2_2, a_10_2_2, a_10_2_2, COMP_CPP) );
00252     INTREPID_TEST_COMMAND( atools.contractDataDataVector<double>(a_10, a_9_2_2, a_10_2_2, COMP_CPP) );
00253     INTREPID_TEST_COMMAND( atools.contractDataDataVector<double>(a_10, a_10_3_2, a_10_2_2, COMP_CPP) );
00254     INTREPID_TEST_COMMAND( atools.contractDataDataVector<double>(a_10, a_10_2_3, a_10_2_2, COMP_CPP) );
00255     INTREPID_TEST_COMMAND( atools.contractDataDataVector<double>(a_2, a_10_2_2, a_10_2_2, COMP_CPP) );
00256     INTREPID_TEST_COMMAND( atools.contractDataDataVector<double>(a_10, a_10_2_2, a_10_2_2, COMP_ENGINE_MAX) );
00257     INTREPID_TEST_COMMAND( atools.contractDataDataVector<double>(a_10, a_10_2_2, a_10_2_2, COMP_CPP) );
00258 
00259     *outStream << "-> contractDataDataTensor:\n";
00260     INTREPID_TEST_COMMAND( atools.contractDataDataTensor<double>(a_2_2, a_2_2, a_2_2, COMP_CPP) );
00261     INTREPID_TEST_COMMAND( atools.contractDataDataTensor<double>(a_2_2, a_10_2_2_2, a_2_2, COMP_CPP) );
00262     INTREPID_TEST_COMMAND( atools.contractDataDataTensor<double>(a_2_2, a_10_2_2_2, a_10_2_2_2, COMP_CPP) );
00263     INTREPID_TEST_COMMAND( atools.contractDataDataTensor<double>(a_10, a_9_2_2_2, a_10_2_2_2, COMP_CPP) );
00264     INTREPID_TEST_COMMAND( atools.contractDataDataTensor<double>(a_10, a_10_2_2_2, a_10_3_2_2, COMP_CPP) );
00265     INTREPID_TEST_COMMAND( atools.contractDataDataTensor<double>(a_10, a_10_2_2_2, a_10_2_3_2, COMP_CPP) );
00266     INTREPID_TEST_COMMAND( atools.contractDataDataTensor<double>(a_10, a_10_2_2_2, a_10_2_2_3, COMP_CPP) );
00267     INTREPID_TEST_COMMAND( atools.contractDataDataTensor<double>(a_2, a_10_2_2_2, a_10_2_2_2, COMP_CPP) );
00268     INTREPID_TEST_COMMAND( atools.contractDataDataTensor<double>(a_10, a_10_2_2_2, a_10_2_2_2, COMP_ENGINE_MAX) );
00269     INTREPID_TEST_COMMAND( atools.contractDataDataTensor<double>(a_10, a_10_2_2_2, a_10_2_2_2, COMP_CPP) );
00270 
00271 #endif
00272 
00273   }
00274   catch (std::logic_error err) {
00275     *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
00276     *outStream << err.what() << '\n';
00277     *outStream << "-------------------------------------------------------------------------------" << "\n\n";
00278     errorFlag = -1000;
00279   };
00280 
00281 #ifdef HAVE_INTREPID_DEBUG
00282   if (Teuchos::TestForException_getThrowNumber() != endThrowNumber)
00283     errorFlag++;
00284 #endif
00285 
00286   *outStream \
00287   << "\n"
00288   << "===============================================================================\n"\
00289   << "| TEST 2: correctness of math operations                                      |\n"\
00290   << "===============================================================================\n";
00291 
00292   outStream->precision(20);
00293 
00294   try {
00295       { // start scope
00296       *outStream << "\n************ Checking contractFieldFieldScalar ************\n";
00297 
00298       int c=5, p=9, l=3, r=7;
00299 
00300       FieldContainer<double> in_c_l_p(c, l, p);
00301       FieldContainer<double> in_c_r_p(c, r, p);
00302       FieldContainer<double> out1_c_l_r(c, l, r);
00303       FieldContainer<double> out2_c_l_r(c, l, r);
00304       double zero = INTREPID_TOL*10000.0;
00305 
00306       // fill with random numbers
00307       for (int i=0; i<in_c_l_p.size(); i++) {
00308         in_c_l_p[i] = Teuchos::ScalarTraits<double>::random();
00309       }
00310       for (int i=0; i<in_c_r_p.size(); i++) {
00311         in_c_r_p[i] = Teuchos::ScalarTraits<double>::random();
00312       }
00313 
00314       art::contractFieldFieldScalar<double>(out1_c_l_r, in_c_l_p, in_c_r_p, COMP_CPP);
00315       art::contractFieldFieldScalar<double>(out2_c_l_r, in_c_l_p, in_c_r_p, COMP_BLAS);
00316       rst::subtract(&out1_c_l_r[0], &out2_c_l_r[0], out2_c_l_r.size());
00317       if (rst::vectorNorm(&out1_c_l_r[0], out1_c_l_r.size(), NORM_ONE) > zero) {
00318         *outStream << "\n\nINCORRECT contractFieldFieldScalar (1): check COMP_CPP vs. COMP_BLAS; "
00319                    << " diff-1norm = " << rst::vectorNorm(&out1_c_l_r[0], out1_c_l_r.size(), NORM_ONE) << "\n\n";
00320         errorFlag++;
00321       }
00322       // with sumInto:
00323       out1_c_l_r.initialize(2.0); out2_c_l_r.initialize(2.0);
00324       art::contractFieldFieldScalar<double>(out1_c_l_r, in_c_l_p, in_c_r_p, COMP_CPP, true);
00325       art::contractFieldFieldScalar<double>(out2_c_l_r, in_c_l_p, in_c_r_p, COMP_BLAS, true);
00326       rst::subtract(&out1_c_l_r[0], &out2_c_l_r[0], out2_c_l_r.size());
00327       if (rst::vectorNorm(&out1_c_l_r[0], out1_c_l_r.size(), NORM_ONE) > zero) {
00328         *outStream << "\n\nINCORRECT contractFieldFieldScalar (1): check COMP_CPP vs. COMP_BLAS; "
00329                    << " diff-1norm = " << rst::vectorNorm(&out1_c_l_r[0], out1_c_l_r.size(), NORM_ONE) << "\n\n";
00330         errorFlag++;
00331       }
00332       } // end scope
00333 
00334       { // start scope
00335       *outStream << "\n************ Checking contractFieldFieldVector ************\n";
00336 
00337       int c=5, p=9, l=3, r=7, d=13;
00338 
00339       FieldContainer<double> in_c_l_p_d(c, l, p, d);
00340       FieldContainer<double> in_c_r_p_d(c, r, p, d);
00341       FieldContainer<double> out1_c_l_r(c, l, r);
00342       FieldContainer<double> out2_c_l_r(c, l, r);
00343       double zero = INTREPID_TOL*10000.0;
00344 
00345       // fill with random numbers
00346       for (int i=0; i<in_c_l_p_d.size(); i++) {
00347         in_c_l_p_d[i] = Teuchos::ScalarTraits<double>::random();
00348       }
00349       for (int i=0; i<in_c_r_p_d.size(); i++) {
00350         in_c_r_p_d[i] = Teuchos::ScalarTraits<double>::random();
00351       }
00352 
00353       art::contractFieldFieldVector<double>(out1_c_l_r, in_c_l_p_d, in_c_r_p_d, COMP_CPP);
00354       art::contractFieldFieldVector<double>(out2_c_l_r, in_c_l_p_d, in_c_r_p_d, COMP_BLAS);
00355 
00356       rst::subtract(&out1_c_l_r[0], &out2_c_l_r[0], out2_c_l_r.size());
00357       if (rst::vectorNorm(&out1_c_l_r[0], out1_c_l_r.size(), NORM_ONE) > zero) {
00358         *outStream << "\n\nINCORRECT contractFieldFieldVector (1): check COMP_CPP vs. COMP_BLAS; "
00359                    << " diff-1norm = " << rst::vectorNorm(&out1_c_l_r[0], out1_c_l_r.size(), NORM_ONE) << "\n\n";
00360         errorFlag++;
00361       }
00362 
00363       // with sumInto:
00364       out1_c_l_r.initialize(2.0); out2_c_l_r.initialize(2.0);
00365 
00366       art::contractFieldFieldVector<double>(out1_c_l_r, in_c_l_p_d, in_c_r_p_d, COMP_CPP, true);
00367       art::contractFieldFieldVector<double>(out2_c_l_r, in_c_l_p_d, in_c_r_p_d, COMP_BLAS, true);
00368 
00369       rst::subtract(&out1_c_l_r[0], &out2_c_l_r[0], out2_c_l_r.size());
00370       if (rst::vectorNorm(&out1_c_l_r[0], out1_c_l_r.size(), NORM_ONE) > zero) {
00371         *outStream << "\n\nINCORRECT contractFieldFieldVector (1): check COMP_CPP vs. COMP_BLAS; "
00372                    << " diff-1norm = " << rst::vectorNorm(&out1_c_l_r[0], out1_c_l_r.size(), NORM_ONE) << "\n\n";
00373         errorFlag++;
00374       }
00375       } // end scope
00376 
00377       { // start scope
00378       *outStream << "\n************ Checking contractFieldFieldTensor ************\n";
00379 
00380       int c=5, p=9, l=3, r=7, d1=13, d2=5;
00381 
00382       FieldContainer<double> in_c_l_p_d_d(c, l, p, d1, d2);
00383       FieldContainer<double> in_c_r_p_d_d(c, r, p, d1, d2);
00384       FieldContainer<double> out1_c_l_r(c, l, r);
00385       FieldContainer<double> out2_c_l_r(c, l, r);
00386       double zero = INTREPID_TOL*10000.0;
00387 
00388       // fill with random numbers
00389       for (int i=0; i<in_c_l_p_d_d.size(); i++) {
00390         in_c_l_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
00391       }
00392       for (int i=0; i<in_c_r_p_d_d.size(); i++) {
00393         in_c_r_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
00394       }
00395 
00396       art::contractFieldFieldTensor<double>(out1_c_l_r, in_c_l_p_d_d, in_c_r_p_d_d, COMP_CPP);
00397       art::contractFieldFieldTensor<double>(out2_c_l_r, in_c_l_p_d_d, in_c_r_p_d_d, COMP_BLAS);
00398 
00399       rst::subtract(&out1_c_l_r[0], &out2_c_l_r[0], out2_c_l_r.size());
00400       if (rst::vectorNorm(&out1_c_l_r[0], out1_c_l_r.size(), NORM_ONE) > zero) {
00401         *outStream << "\n\nINCORRECT contractFieldFieldTensor (1): check COMP_CPP vs. COMP_BLAS; "
00402                    << " diff-1norm = " << rst::vectorNorm(&out1_c_l_r[0], out1_c_l_r.size(), NORM_ONE) << "\n\n";
00403         errorFlag++;
00404       }
00405 
00406       // with sumInto:
00407       out1_c_l_r.initialize(2.0); out2_c_l_r.initialize(2.0);
00408 
00409       art::contractFieldFieldTensor<double>(out1_c_l_r, in_c_l_p_d_d, in_c_r_p_d_d, COMP_CPP, true);
00410       art::contractFieldFieldTensor<double>(out2_c_l_r, in_c_l_p_d_d, in_c_r_p_d_d, COMP_BLAS, true);
00411 
00412       rst::subtract(&out1_c_l_r[0], &out2_c_l_r[0], out2_c_l_r.size());
00413       if (rst::vectorNorm(&out1_c_l_r[0], out1_c_l_r.size(), NORM_ONE) > zero) {
00414         *outStream << "\n\nINCORRECT contractFieldFieldTensor (1): check COMP_CPP vs. COMP_BLAS; "
00415                    << " diff-1norm = " << rst::vectorNorm(&out1_c_l_r[0], out1_c_l_r.size(), NORM_ONE) << "\n\n";
00416         errorFlag++;
00417       }
00418       } // end scope
00419 
00420       { // start scope
00421       *outStream << "\n************ Checking contractDataFieldScalar ************\n";
00422 
00423       int c=5, p=9, l=7;
00424 
00425       FieldContainer<double> in_c_l_p(c, l, p);
00426       FieldContainer<double> data_c_p(c, p);
00427       FieldContainer<double> data_c_1(c, 1);
00428       FieldContainer<double> out1_c_l(c, l);
00429       FieldContainer<double> out2_c_l(c, l);
00430       double zero = INTREPID_TOL*10000.0;
00431 
00432       // fill with random numbers
00433       for (int i=0; i<in_c_l_p.size(); i++) {
00434         in_c_l_p[i] = Teuchos::ScalarTraits<double>::random();
00435       }
00436       for (int i=0; i<data_c_p.size(); i++) {
00437         data_c_p[i] = Teuchos::ScalarTraits<double>::random();
00438       }
00439       for (int i=0; i<data_c_1.size(); i++) {
00440         data_c_1[i] = Teuchos::ScalarTraits<double>::random();
00441       }
00442 
00443       // nonconstant data
00444       art::contractDataFieldScalar<double>(out1_c_l, data_c_p, in_c_l_p, COMP_CPP);
00445       art::contractDataFieldScalar<double>(out2_c_l, data_c_p, in_c_l_p, COMP_BLAS);
00446       rst::subtract(&out1_c_l[0], &out2_c_l[0], out2_c_l.size());
00447       if (rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) > zero) {
00448         *outStream << "\n\nINCORRECT contractDataFieldScalar (1): check COMP_CPP vs. COMP_BLAS; "
00449                    << " diff-1norm = " << rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) << "\n\n";
00450         errorFlag++;
00451       }
00452       // constant data
00453       art::contractDataFieldScalar<double>(out1_c_l, data_c_1, in_c_l_p, COMP_CPP);
00454       art::contractDataFieldScalar<double>(out2_c_l, data_c_1, in_c_l_p, COMP_BLAS);
00455       rst::subtract(&out1_c_l[0], &out2_c_l[0], out2_c_l.size());
00456       if (rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) > zero) {
00457         *outStream << "\n\nINCORRECT contractDataFieldScalar (2): check COMP_CPP vs. COMP_BLAS; "
00458                    << " diff-1norm = " << rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) << "\n\n";
00459         errorFlag++;
00460       }
00461       // nonconstant data with sumInto
00462       out1_c_l.initialize(2.0); out2_c_l.initialize(2.0);
00463       art::contractDataFieldScalar<double>(out1_c_l, data_c_p, in_c_l_p, COMP_CPP, true);
00464       art::contractDataFieldScalar<double>(out2_c_l, data_c_p, in_c_l_p, COMP_BLAS, true);
00465       rst::subtract(&out1_c_l[0], &out2_c_l[0], out2_c_l.size());
00466       if (rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) > zero) {
00467         *outStream << "\n\nINCORRECT contractDataFieldScalar (1): check COMP_CPP vs. COMP_BLAS; "
00468                    << " diff-1norm = " << rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) << "\n\n";
00469         errorFlag++;
00470       }
00471       } // end scope
00472 
00473       { // start scope
00474       *outStream << "\n************ Checking contractDataFieldVector ************\n";
00475 
00476       int c=5, p=9, l=7, d=3;
00477 
00478       FieldContainer<double> in_c_l_p_d(c, l, p, d);
00479       FieldContainer<double> data_c_p_d(c, p, d);
00480       FieldContainer<double> data_c_1_d(c, 1, d);
00481       FieldContainer<double> out1_c_l(c, l);
00482       FieldContainer<double> out2_c_l(c, l);
00483       double zero = INTREPID_TOL*10000.0;
00484 
00485       // fill with random numbers
00486       for (int i=0; i<in_c_l_p_d.size(); i++) {
00487         in_c_l_p_d[i] = Teuchos::ScalarTraits<double>::random();
00488       }
00489       for (int i=0; i<data_c_p_d.size(); i++) {
00490         data_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
00491       }
00492       for (int i=0; i<data_c_1_d.size(); i++) {
00493         data_c_1_d[i] = Teuchos::ScalarTraits<double>::random();
00494       }
00495 
00496       // nonconstant data
00497       art::contractDataFieldVector<double>(out1_c_l, data_c_p_d, in_c_l_p_d, COMP_CPP);
00498       art::contractDataFieldVector<double>(out2_c_l, data_c_p_d, in_c_l_p_d, COMP_BLAS);
00499       rst::subtract(&out1_c_l[0], &out2_c_l[0], out2_c_l.size());
00500       if (rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) > zero) {
00501         *outStream << "\n\nINCORRECT contractDataFieldVector (1): check COMP_CPP vs. COMP_BLAS; "
00502                    << " diff-1norm = " << rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) << "\n\n";
00503         errorFlag++;
00504       }
00505       // constant data
00506       art::contractDataFieldVector<double>(out1_c_l, data_c_1_d, in_c_l_p_d, COMP_CPP);
00507       art::contractDataFieldVector<double>(out2_c_l, data_c_1_d, in_c_l_p_d, COMP_BLAS);
00508       rst::subtract(&out1_c_l[0], &out2_c_l[0], out2_c_l.size());
00509       if (rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) > zero) {
00510         *outStream << "\n\nINCORRECT contractDataFieldVector (2): check COMP_CPP vs. COMP_BLAS; "
00511                    << " diff-1norm = " << rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) << "\n\n";
00512         errorFlag++;
00513       }
00514       // nonconstant data with sumInto
00515       out1_c_l.initialize(2.0); out2_c_l.initialize(2.0);
00516       art::contractDataFieldVector<double>(out1_c_l, data_c_p_d, in_c_l_p_d, COMP_CPP, true);
00517       art::contractDataFieldVector<double>(out2_c_l, data_c_p_d, in_c_l_p_d, COMP_BLAS, true);
00518       rst::subtract(&out1_c_l[0], &out2_c_l[0], out2_c_l.size());
00519       if (rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) > zero) {
00520         *outStream << "\n\nINCORRECT contractDataFieldVector (1): check COMP_CPP vs. COMP_BLAS; "
00521                    << " diff-1norm = " << rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) << "\n\n";
00522         errorFlag++;
00523       }
00524       } // end scope
00525 
00526       { // start scope
00527       *outStream << "\n************ Checking contractDataFieldTensor ************\n";
00528 
00529       int c=5, p=9, l=7, d1=3, d2=13;
00530 
00531       FieldContainer<double> in_c_l_p_d_d(c, l, p, d1, d2);
00532       FieldContainer<double> data_c_p_d_d(c, p, d1, d2);
00533       FieldContainer<double> data_c_1_d_d(c, 1, d1, d2);
00534       FieldContainer<double> out1_c_l(c, l);
00535       FieldContainer<double> out2_c_l(c, l);
00536       double zero = INTREPID_TOL*10000.0;
00537 
00538       // fill with random numbers
00539       for (int i=0; i<in_c_l_p_d_d.size(); i++) {
00540         in_c_l_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
00541       }
00542       for (int i=0; i<data_c_p_d_d.size(); i++) {
00543         data_c_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
00544       }
00545       for (int i=0; i<data_c_1_d_d.size(); i++) {
00546         data_c_1_d_d[i] = Teuchos::ScalarTraits<double>::random();
00547       }
00548 
00549       // nonconstant data
00550       art::contractDataFieldTensor<double>(out1_c_l, data_c_p_d_d, in_c_l_p_d_d, COMP_CPP);
00551       art::contractDataFieldTensor<double>(out2_c_l, data_c_p_d_d, in_c_l_p_d_d, COMP_BLAS);
00552       rst::subtract(&out1_c_l[0], &out2_c_l[0], out2_c_l.size());
00553       if (rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) > zero) {
00554         *outStream << "\n\nINCORRECT contractDataFieldTensor (1): check COMP_CPP vs. COMP_BLAS; "
00555                    << " diff-1norm = " << rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) << "\n\n";
00556         errorFlag++;
00557       }
00558       // constant data
00559       art::contractDataFieldTensor<double>(out1_c_l, data_c_1_d_d, in_c_l_p_d_d, COMP_CPP);
00560       art::contractDataFieldTensor<double>(out2_c_l, data_c_1_d_d, in_c_l_p_d_d, COMP_BLAS);
00561       rst::subtract(&out1_c_l[0], &out2_c_l[0], out2_c_l.size());
00562       if (rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) > zero) {
00563         *outStream << "\n\nINCORRECT contractDataFieldTensor (1): check COMP_CPP vs. COMP_BLAS; "
00564                    << " diff-1norm = " << rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) << "\n\n";
00565         errorFlag++;
00566       }
00567       // nonconstant data with sumInto
00568       out1_c_l.initialize(2.0); out2_c_l.initialize(2.0);
00569       art::contractDataFieldTensor<double>(out1_c_l, data_c_p_d_d, in_c_l_p_d_d, COMP_CPP, true);
00570       art::contractDataFieldTensor<double>(out2_c_l, data_c_p_d_d, in_c_l_p_d_d, COMP_BLAS, true);
00571       rst::subtract(&out1_c_l[0], &out2_c_l[0], out2_c_l.size());
00572       if (rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) > zero) {
00573         *outStream << "\n\nINCORRECT contractDataFieldTensor (1): check COMP_CPP vs. COMP_BLAS; "
00574                    << " diff-1norm = " << rst::vectorNorm(&out1_c_l[0], out1_c_l.size(), NORM_ONE) << "\n\n";
00575         errorFlag++;
00576       }
00577       } // end scope
00578 
00579       { // start scope
00580       *outStream << "\n************ Checking contractDataDataScalar ************\n";
00581 
00582       int c=5, p=9;
00583 
00584       FieldContainer<double> inl_c_p(c, p);
00585       FieldContainer<double> inr_c_p(c, p);
00586       FieldContainer<double> out1_c(c);
00587       FieldContainer<double> out2_c(c);
00588       double zero = INTREPID_TOL*10000.0;
00589 
00590       // fill with random numbers
00591       for (int i=0; i<inl_c_p.size(); i++) {
00592         inl_c_p[i] = Teuchos::ScalarTraits<double>::random();
00593       }
00594       for (int i=0; i<inr_c_p.size(); i++) {
00595         inr_c_p[i] = Teuchos::ScalarTraits<double>::random();
00596       }
00597 
00598       art::contractDataDataScalar<double>(out1_c, inl_c_p, inr_c_p, COMP_CPP);
00599       art::contractDataDataScalar<double>(out2_c, inl_c_p, inr_c_p, COMP_BLAS);
00600       rst::subtract(&out1_c[0], &out2_c[0], out2_c.size());
00601       if (rst::vectorNorm(&out1_c[0], out1_c.size(), NORM_ONE) > zero) {
00602         *outStream << "\n\nINCORRECT contractDataDataScalar (1): check COMP_CPP vs. COMP_BLAS; "
00603                    << " diff-1norm = " << rst::vectorNorm(&out1_c[0], out1_c.size(), NORM_ONE) << "\n\n";
00604         errorFlag++;
00605       }
00606       // with sumInto:
00607       out1_c.initialize(2.0); out2_c.initialize(2.0);
00608       art::contractDataDataScalar<double>(out1_c, inl_c_p, inr_c_p, COMP_CPP, true);
00609       art::contractDataDataScalar<double>(out2_c, inl_c_p, inr_c_p, COMP_BLAS, true);
00610       rst::subtract(&out1_c[0], &out2_c[0], out2_c.size());
00611       if (rst::vectorNorm(&out1_c[0], out1_c.size(), NORM_ONE) > zero) {
00612         *outStream << "\n\nINCORRECT contractDataDataScalar (1): check COMP_CPP vs. COMP_BLAS; "
00613                    << " diff-1norm = " << rst::vectorNorm(&out1_c[0], out1_c.size(), NORM_ONE) << "\n\n";
00614         errorFlag++;
00615       }
00616       } // end scope
00617 
00618       { // start scope
00619       *outStream << "\n************ Checking contractDataDataVector ************\n";
00620 
00621       int c=5, p=9, d=13;
00622 
00623       FieldContainer<double> inl_c_p_d(c, p, d);
00624       FieldContainer<double> inr_c_p_d(c, p, d);
00625       FieldContainer<double> out1_c(c);
00626       FieldContainer<double> out2_c(c);
00627       double zero = INTREPID_TOL*10000.0;
00628 
00629       // fill with random numbers
00630       for (int i=0; i<inl_c_p_d.size(); i++) {
00631         inl_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
00632       }
00633       for (int i=0; i<inr_c_p_d.size(); i++) {
00634         inr_c_p_d[i] = Teuchos::ScalarTraits<double>::random();
00635       }
00636 
00637       art::contractDataDataVector<double>(out1_c, inl_c_p_d, inr_c_p_d, COMP_CPP);
00638       art::contractDataDataVector<double>(out2_c, inl_c_p_d, inr_c_p_d, COMP_BLAS);
00639 
00640       rst::subtract(&out1_c[0], &out2_c[0], out2_c.size());
00641       if (rst::vectorNorm(&out1_c[0], out1_c.size(), NORM_ONE) > zero) {
00642         *outStream << "\n\nINCORRECT contractDataDataVector (1): check COMP_CPP vs. COMP_BLAS; "
00643                    << " diff-1norm = " << rst::vectorNorm(&out1_c[0], out1_c.size(), NORM_ONE) << "\n\n";
00644         errorFlag++;
00645       }
00646 
00647       // with sumInto:
00648       out1_c.initialize(2.0); out2_c.initialize(2.0);
00649 
00650       art::contractDataDataVector<double>(out1_c, inl_c_p_d, inr_c_p_d, COMP_CPP, true);
00651       art::contractDataDataVector<double>(out2_c, inl_c_p_d, inr_c_p_d, COMP_BLAS, true);
00652 
00653       rst::subtract(&out1_c[0], &out2_c[0], out2_c.size());
00654       if (rst::vectorNorm(&out1_c[0], out1_c.size(), NORM_ONE) > zero) {
00655         *outStream << "\n\nINCORRECT contractDataDataVector (1): check COMP_CPP vs. COMP_BLAS; "
00656                    << " diff-1norm = " << rst::vectorNorm(&out1_c[0], out1_c.size(), NORM_ONE) << "\n\n";
00657         errorFlag++;
00658       }
00659       } // end scope
00660 
00661       { // start scope
00662       *outStream << "\n************ Checking contractDataDataTensor ************\n";
00663 
00664       int c=5, p=9, d1=13, d2=5;
00665 
00666       FieldContainer<double> inl_c_p_d_d(c, p, d1, d2);
00667       FieldContainer<double> inr_c_p_d_d(c, p, d1, d2);
00668       FieldContainer<double> out1_c(c);
00669       FieldContainer<double> out2_c(c);
00670       double zero = INTREPID_TOL*10000.0;
00671 
00672       // fill with random numbers
00673       for (int i=0; i<inl_c_p_d_d.size(); i++) {
00674         inl_c_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
00675       }
00676       for (int i=0; i<inr_c_p_d_d.size(); i++) {
00677         inr_c_p_d_d[i] = Teuchos::ScalarTraits<double>::random();
00678       }
00679 
00680       art::contractDataDataTensor<double>(out1_c, inl_c_p_d_d, inr_c_p_d_d, COMP_CPP);
00681       art::contractDataDataTensor<double>(out2_c, inl_c_p_d_d, inr_c_p_d_d, COMP_BLAS);
00682 
00683       rst::subtract(&out1_c[0], &out2_c[0], out2_c.size());
00684       if (rst::vectorNorm(&out1_c[0], out1_c.size(), NORM_ONE) > zero) {
00685         *outStream << "\n\nINCORRECT contractDataDataTensor (1): check COMP_CPP vs. COMP_BLAS; "
00686                    << " diff-1norm = " << rst::vectorNorm(&out1_c[0], out1_c.size(), NORM_ONE) << "\n\n";
00687         errorFlag++;
00688       }
00689 
00690       // with sumInto:
00691       out1_c.initialize(2.0); out2_c.initialize(2.0);
00692 
00693       art::contractDataDataTensor<double>(out1_c, inl_c_p_d_d, inr_c_p_d_d, COMP_CPP, true);
00694       art::contractDataDataTensor<double>(out2_c, inl_c_p_d_d, inr_c_p_d_d, COMP_BLAS, true);
00695 
00696       rst::subtract(&out1_c[0], &out2_c[0], out2_c.size());
00697       if (rst::vectorNorm(&out1_c[0], out1_c.size(), NORM_ONE) > zero) {
00698         *outStream << "\n\nINCORRECT contractDataDataTensor (1): check COMP_CPP vs. COMP_BLAS; "
00699                    << " diff-1norm = " << rst::vectorNorm(&out1_c[0], out1_c.size(), NORM_ONE) << "\n\n";
00700         errorFlag++;
00701       }
00702       } // end scope
00703 
00704       /******************************************/
00705       *outStream << "\n";
00706   }
00707   catch (std::logic_error err) {
00708     *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
00709     *outStream << err.what() << '\n';
00710     *outStream << "-------------------------------------------------------------------------------" << "\n\n";
00711     errorFlag = -1000;
00712   };
00713 
00714 
00715   if (errorFlag != 0)
00716     std::cout << "End Result: TEST FAILED\n";
00717   else
00718     std::cout << "End Result: TEST PASSED\n";
00719 
00720   // reset format state of std::cout
00721   std::cout.copyfmt(oldFormatState);
00722 
00723   return errorFlag;
00724 }