Intrepid
http://trilinos.sandia.gov/packages/docs/r10.10/packages/intrepid/test/Shared/ArrayTools/test_03.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: dot multiply                                       |\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 + 36;
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);
00129     FieldContainer<double> a_2_2(2, 2);
00130     FieldContainer<double> a_3_2(3, 2);
00131     FieldContainer<double> a_2_2_2(2, 2, 2);
00132     FieldContainer<double> a_2_2_2_2(2, 2, 2, 2);
00133     FieldContainer<double> a_10_1(10, 1);
00134     FieldContainer<double> a_10_2(10, 2);
00135     FieldContainer<double> a_10_3(10, 3);
00136     FieldContainer<double> a_10_1_2(10, 1, 2);
00137     FieldContainer<double> a_10_2_2(10, 2, 2);
00138     FieldContainer<double> a_10_2_2_2(10, 2, 2, 2);
00139     FieldContainer<double> a_10_2_2_2_2(10, 2, 2, 2, 2);
00140     FieldContainer<double> a_9_2_2(9, 2, 2);
00141     FieldContainer<double> a_10_3_2(10, 3, 2);
00142     FieldContainer<double> a_10_2_3(10, 2, 3);
00143     FieldContainer<double> a_10_2_2_3(10, 2, 2, 3);
00144     FieldContainer<double> a_10_1_2_2(10, 1, 2, 2);
00145 
00146     *outStream << "-> dotMultiplyDataField:\n";
00147     INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_2_2, a_2, a_2_2) );
00148     INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_2_2, a_10_2_2, a_10_2_2_2_2) );
00149     INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2, a_10_2_2, a_10_2_2_2) );
00150     INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_9_2_2, a_10_2_2_2) );
00151     INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_10_3_2, a_10_2_2_2) );
00152     INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_10_2_3, a_10_2_2_2) );
00153     INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_10_2_2_3, a_10_2_2_2_2) );
00154     INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_2_2_2, a_10_2_2_2, a_10_2_2_2_2) );
00155     INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_3_2, a_10_2_2_2, a_10_2_2_2_2) );
00156     INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_3, a_10_2_2_2, a_10_2_2_2_2) );
00157     INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_10_2_2_2, a_10_2_2_2_2) );
00158     INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_10_1_2_2, a_10_2_2_2_2) );
00159     //
00160     INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_2_2, a_2, a_2) );
00161     INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_2_2, a_10_2_2, a_10_2) );
00162     INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2, a_10_2_2, a_10_2_2) );
00163     INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_10_2_2, a_10_3_2) );
00164     INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_10_2_2, a_10_2_2) );
00165     INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_3, a_10_2_2, a_2_2_2) );
00166     INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_9_2_2, a_2_2_2) );
00167     INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_10_2_3, a_2_2_2) );
00168     INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_10_2_2_3, a_2_2_2_2) );
00169     INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_10_2_2_2, a_2_2_2_2) );
00170     INTREPID_TEST_COMMAND( atools.dotMultiplyDataField<double>(a_10_2_2, a_10_1, a_2_2) );
00171 
00172     *outStream << "-> dotMultiplyDataData:\n";
00173     INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_2_2, a_2, a_2_2) );
00174     INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_2_2, a_10_2_2, a_10_2_2_2) );
00175     INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2_2, a_10_2_2, a_10_2_2) );
00176     INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_9_2_2, a_10_2_2) );
00177     INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_3_2, a_10_2_2) );
00178     INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_2_3, a_10_2_2) );
00179     INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_2_2_2, a_10_2_2_3) );
00180     INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_9_2_2, a_9_2_2) );
00181     INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_3_2, a_10_3_2) );
00182     INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_2, a_10_2) );
00183     INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_2_2, a_10_2_2) );
00184     INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_2_2_2, a_10_2_2_2) );
00185     INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_1, a_10_2) );
00186     INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_1_2, a_10_2_2) );
00187     INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_1_2_2, a_10_2_2_2) );
00188     //
00189     INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_2_2, a_10_2_2_2_2, a_10_2) );
00190     INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_2_2, a_10_2_2, a_2) );
00191     INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2_2, a_10_2_2, a_10_2) );
00192     INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_2_2, a_10_2) );
00193     INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_3, a_10_2_2, a_2_2) );
00194     INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_9_2_2, a_2_2) );
00195     INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_2_3, a_2_2) );
00196     INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_2_2_3, a_2_2_2) );
00197     INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_2, a_2) );
00198     INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_2_2, a_2_2) );
00199     INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_2_2_2, a_2_2_2) );
00200     INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_1, a_2) );
00201     INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_1_2, a_2_2) );
00202     INTREPID_TEST_COMMAND( atools.dotMultiplyDataData<double>(a_10_2, a_10_1_2_2, a_2_2_2) );
00203 #endif
00204 
00205   }
00206   catch (std::logic_error err) {
00207     *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
00208     *outStream << err.what() << '\n';
00209     *outStream << "-------------------------------------------------------------------------------" << "\n\n";
00210     errorFlag = -1000;
00211   };
00212 
00213 #ifdef HAVE_INTREPID_DEBUG
00214   if (Teuchos::TestForException_getThrowNumber() != endThrowNumber)
00215     errorFlag++;
00216 #endif
00217 
00218   *outStream \
00219   << "\n"
00220   << "===============================================================================\n"\
00221   << "| TEST 2: correctness of math operations                                      |\n"\
00222   << "===============================================================================\n";
00223 
00224   outStream->precision(20);
00225 
00226   try {
00227       { // start scope
00228       *outStream << "\n************ Checking dotMultiplyDataField, (branch 1) ************\n";
00229 
00230       int c=5, p=9, f=7, d1=6, d2=14;
00231 
00232       FieldContainer<double> in_c_f_p(c, f, p);
00233       FieldContainer<double> in_c_f_p_d(c, f, p, d1);
00234       FieldContainer<double> in_c_f_p_d_d(c, f, p, d1, d2);
00235       FieldContainer<double> data_c_p(c, p);
00236       FieldContainer<double> data_c_p_d(c, p, d1);
00237       FieldContainer<double> data_c_p_d_d(c, p, d1, d2);
00238       FieldContainer<double> data_c_1(c, 1);
00239       FieldContainer<double> data_c_1_d(c, 1, d1);
00240       FieldContainer<double> data_c_1_d_d(c, 1, d1, d2);
00241       FieldContainer<double> out_c_f_p(c, f, p);
00242       FieldContainer<double> outSM_c_f_p(c, f, p);
00243       FieldContainer<double> outDM_c_f_p(c, f, p);
00244       double zero = INTREPID_TOL*10000.0;
00245 
00246       // fill with random numbers
00247       for (int i=0; i<in_c_f_p.size(); i++) {
00248         in_c_f_p[i] = Teuchos::ScalarTraits<double>::random();
00249       }
00250       // fill with alternating 1's and -1's
00251       for (int i=0; i<in_c_f_p_d.size(); i++) {
00252         in_c_f_p_d[i] = std::pow((double)-1.0, (int)i);
00253       }
00254       // fill with 1's
00255       for (int i=0; i<in_c_f_p_d_d.size(); i++) {
00256         in_c_f_p_d_d[i] = 1.0;
00257       }
00258       // fill with random numbers
00259       for (int i=0; i<data_c_p.size(); i++) {
00260         data_c_p[i] = Teuchos::ScalarTraits<double>::random();
00261       }
00262       // fill with 1's
00263       for (int i=0; i<data_c_1.size(); i++) {
00264         data_c_1[i] = Teuchos::ScalarTraits<double>::random();
00265       }
00266       // fill with 1's
00267       for (int i=0; i<data_c_p_d.size(); i++) {
00268         data_c_p_d[i] = 1.0;
00269       }
00270       // fill with 1's
00271       for (int i=0; i<data_c_1_d.size(); i++) {
00272         data_c_1_d[i] = 1.0;
00273       }
00274       // fill with 1's
00275       for (int i=0; i<data_c_p_d_d.size(); i++) {
00276         data_c_p_d_d[i] = 1.0;
00277       }
00278       // fill with 1's
00279       for (int i=0; i<data_c_1_d_d.size(); i++) {
00280         data_c_1_d_d[i] = 1.0;
00281       }
00282 
00283       art::scalarMultiplyDataField<double>(outSM_c_f_p, data_c_p, in_c_f_p);
00284       art::dotMultiplyDataField<double>(outDM_c_f_p, data_c_p, in_c_f_p);
00285       rst::subtract(out_c_f_p, outSM_c_f_p, outDM_c_f_p);
00286       if (rst::vectorNorm(&out_c_f_p[0], out_c_f_p.size(), NORM_ONE) > zero) {
00287         *outStream << "\n\nINCORRECT dotMultiplyDataField (1): check dot multiply for scalars vs. scalar multiply\n\n";
00288         errorFlag = -1000;
00289       }
00290       art::dotMultiplyDataField<double>(out_c_f_p, data_c_p_d, in_c_f_p_d);
00291       if (rst::vectorNorm(&out_c_f_p[0], out_c_f_p.size(), NORM_ONE) > zero) {
00292         *outStream << "\n\nINCORRECT dotMultiplyDataField (2): check dot multiply of orthogonal vectors\n\n";
00293         errorFlag = -1000;
00294       }
00295       art::dotMultiplyDataField<double>(out_c_f_p, data_c_p_d_d, in_c_f_p_d_d);
00296       if ((rst::vectorNorm(&out_c_f_p[0], out_c_f_p.size(), NORM_INF) - d1*d2) > zero) {
00297         *outStream << "\n\nINCORRECT dotMultiplyDataField (3): check dot multiply for tensors of 1s\n\n";
00298         errorFlag = -1000;
00299       }
00300       art::scalarMultiplyDataField<double>(outSM_c_f_p, data_c_1, in_c_f_p);
00301       art::dotMultiplyDataField<double>(outDM_c_f_p, data_c_1, in_c_f_p);
00302       rst::subtract(out_c_f_p, outSM_c_f_p, outDM_c_f_p);
00303       if (rst::vectorNorm(&out_c_f_p[0], out_c_f_p.size(), NORM_ONE) > zero) {
00304         *outStream << "\n\nINCORRECT dotMultiplyDataField (4): check dot multiply for scalars vs. scalar multiply\n\n";
00305         errorFlag = -1000;
00306       }
00307       art::dotMultiplyDataField<double>(out_c_f_p, data_c_1_d, in_c_f_p_d);
00308       if (rst::vectorNorm(&out_c_f_p[0], out_c_f_p.size(), NORM_ONE) > zero) {
00309         *outStream << "\n\nINCORRECT dotMultiplyDataField (5): check dot multiply of orthogonal vectors\n\n";
00310         errorFlag = -1000;
00311       }
00312       art::dotMultiplyDataField<double>(out_c_f_p, data_c_1_d_d, in_c_f_p_d_d);
00313       if ((rst::vectorNorm(&out_c_f_p[0], out_c_f_p.size(), NORM_INF) - d1*d2) > zero) {
00314         *outStream << "\n\nINCORRECT dotMultiplyDataField (6): check dot multiply for tensors of 1s\n\n";
00315         errorFlag = -1000;
00316       }
00317       } // end scope
00318 
00319 
00320       { // start scope
00321       *outStream << "\n************ Checking dotMultiplyDataField, (branch 2) ************\n";
00322 
00323       int c=5, p=9, f=7, d1=6, d2=14;
00324 
00325       FieldContainer<double> in_f_p(f, p);
00326       FieldContainer<double> in_f_p_d(f, p, d1);
00327       FieldContainer<double> in_f_p_d_d(f, p, d1, d2);
00328       FieldContainer<double> data_c_p(c, p);
00329       FieldContainer<double> data_c_p_d(c, p, d1);
00330       FieldContainer<double> data_c_p_d_d(c, p, d1, d2);
00331       FieldContainer<double> data_c_1(c, 1);
00332       FieldContainer<double> data_c_1_d(c, 1, d1);
00333       FieldContainer<double> data_c_1_d_d(c, 1, d1, d2);
00334       FieldContainer<double> out_c_f_p(c, f, p);
00335       FieldContainer<double> outSM_c_f_p(c, f, p);
00336       FieldContainer<double> outDM_c_f_p(c, f, p);
00337       double zero = INTREPID_TOL*10000.0;
00338 
00339       // fill with random numbers
00340       for (int i=0; i<in_f_p.size(); i++) {
00341         in_f_p[i] = Teuchos::ScalarTraits<double>::random();
00342       }
00343       // fill with alternating 1's and -1's
00344       for (int i=0; i<in_f_p_d.size(); i++) {
00345         in_f_p_d[i] = std::pow((double)-1.0, (int)i);
00346       }
00347       // fill with 1's
00348       for (int i=0; i<in_f_p_d_d.size(); i++) {
00349         in_f_p_d_d[i] = 1.0;
00350       }
00351       // fill with random numbers
00352       for (int i=0; i<data_c_p.size(); i++) {
00353         data_c_p[i] = Teuchos::ScalarTraits<double>::random();
00354       }
00355       // fill with 1's
00356       for (int i=0; i<data_c_1.size(); i++) {
00357         data_c_1[i] = Teuchos::ScalarTraits<double>::random();
00358       }
00359       // fill with 1's
00360       for (int i=0; i<data_c_p_d.size(); i++) {
00361         data_c_p_d[i] = 1.0;
00362       }
00363       // fill with 1's
00364       for (int i=0; i<data_c_1_d.size(); i++) {
00365         data_c_1_d[i] = 1.0;
00366       }
00367       // fill with 1's
00368       for (int i=0; i<data_c_p_d_d.size(); i++) {
00369         data_c_p_d_d[i] = 1.0;
00370       }
00371       // fill with 1's
00372       for (int i=0; i<data_c_1_d_d.size(); i++) {
00373         data_c_1_d_d[i] = 1.0;
00374       }
00375 
00376       art::scalarMultiplyDataField<double>(outSM_c_f_p, data_c_p, in_f_p);
00377       art::dotMultiplyDataField<double>(outDM_c_f_p, data_c_p, in_f_p);
00378       rst::subtract(out_c_f_p, outSM_c_f_p, outDM_c_f_p);
00379       if (rst::vectorNorm(&out_c_f_p[0], out_c_f_p.size(), NORM_ONE) > zero) {
00380         *outStream << "\n\nINCORRECT dotMultiplyDataField (7): check dot multiply for scalars vs. scalar multiply\n\n";
00381         errorFlag = -1000;
00382       }
00383       art::dotMultiplyDataField<double>(out_c_f_p, data_c_p_d, in_f_p_d);
00384       if (rst::vectorNorm(&out_c_f_p[0], out_c_f_p.size(), NORM_ONE) > zero) {
00385         *outStream << "\n\nINCORRECT dotMultiplyDataField (8): check dot multiply of orthogonal vectors\n\n";
00386         errorFlag = -1000;
00387       }
00388       art::dotMultiplyDataField<double>(out_c_f_p, data_c_p_d_d, in_f_p_d_d);
00389       if ((rst::vectorNorm(&out_c_f_p[0], out_c_f_p.size(), NORM_INF) - d1*d2) > zero) {
00390         *outStream << "\n\nINCORRECT dotMultiplyDataField (9): check dot multiply for tensors of 1s\n\n";
00391         errorFlag = -1000;
00392       }
00393       art::scalarMultiplyDataField<double>(outSM_c_f_p, data_c_1, in_f_p);
00394       art::dotMultiplyDataField<double>(outDM_c_f_p, data_c_1, in_f_p);
00395       rst::subtract(out_c_f_p, outSM_c_f_p, outDM_c_f_p);
00396       if (rst::vectorNorm(&out_c_f_p[0], out_c_f_p.size(), NORM_ONE) > zero) {
00397         *outStream << "\n\nINCORRECT dotMultiplyDataField (10): check dot multiply for scalars vs. scalar multiply\n\n";
00398         errorFlag = -1000;
00399       }
00400       art::dotMultiplyDataField<double>(out_c_f_p, data_c_1_d, in_f_p_d);
00401       if (rst::vectorNorm(&out_c_f_p[0], out_c_f_p.size(), NORM_ONE) > zero) {
00402         *outStream << "\n\nINCORRECT dotMultiplyDataField (11): check dot multiply of orthogonal vectors\n\n";
00403         errorFlag = -1000;
00404       }
00405       art::dotMultiplyDataField<double>(out_c_f_p, data_c_1_d_d, in_f_p_d_d);
00406       if ((rst::vectorNorm(&out_c_f_p[0], out_c_f_p.size(), NORM_INF) - d1*d2) > zero) {
00407         *outStream << "\n\nINCORRECT dotMultiplyDataField (12): check dot multiply for tensors of 1s\n\n";
00408         errorFlag = -1000;
00409       }
00410       } // end scope
00411 
00412 
00413 
00414 
00415       { // start scope
00416       *outStream << "\n************ Checking dotMultiplyDataData, (branch 1) ************\n";
00417 
00418       int c=5, p=9, d1=6, d2=14;
00419 
00420       FieldContainer<double> in_c_p(c, p);
00421       FieldContainer<double> in_c_p_d(c, p, d1);
00422       FieldContainer<double> in_c_p_d_d(c, p, d1, d2);
00423       FieldContainer<double> data_c_p(c, p);
00424       FieldContainer<double> data_c_p_d(c, p, d1);
00425       FieldContainer<double> data_c_p_d_d(c, p, d1, d2);
00426       FieldContainer<double> data_c_1(c, 1);
00427       FieldContainer<double> data_c_1_d(c, 1, d1);
00428       FieldContainer<double> data_c_1_d_d(c, 1, d1, d2);
00429       FieldContainer<double> out_c_p(c, p);
00430       FieldContainer<double> outSM_c_p(c, p);
00431       FieldContainer<double> outDM_c_p(c, p);
00432       double zero = INTREPID_TOL*10000.0;
00433 
00434       // fill with random numbers
00435       for (int i=0; i<in_c_p.size(); i++) {
00436         in_c_p[i] = Teuchos::ScalarTraits<double>::random();
00437       }
00438       // fill with alternating 1's and -1's
00439       for (int i=0; i<in_c_p_d.size(); i++) {
00440         in_c_p_d[i] = std::pow((double)-1.0, (int)i);
00441       }
00442       // fill with 1's
00443       for (int i=0; i<in_c_p_d_d.size(); i++) {
00444         in_c_p_d_d[i] = 1.0;
00445       }
00446       // fill with random numbers
00447       for (int i=0; i<data_c_p.size(); i++) {
00448         data_c_p[i] = Teuchos::ScalarTraits<double>::random();
00449       }
00450       // fill with 1's
00451       for (int i=0; i<data_c_1.size(); i++) {
00452         data_c_1[i] = Teuchos::ScalarTraits<double>::random();
00453       }
00454       // fill with 1's
00455       for (int i=0; i<data_c_p_d.size(); i++) {
00456         data_c_p_d[i] = 1.0;
00457       }
00458       // fill with 1's
00459       for (int i=0; i<data_c_1_d.size(); i++) {
00460         data_c_1_d[i] = 1.0;
00461       }
00462       // fill with 1's
00463       for (int i=0; i<data_c_p_d_d.size(); i++) {
00464         data_c_p_d_d[i] = 1.0;
00465       }
00466       // fill with 1's
00467       for (int i=0; i<data_c_1_d_d.size(); i++) {
00468         data_c_1_d_d[i] = 1.0;
00469       }
00470 
00471       art::scalarMultiplyDataData<double>(outSM_c_p, data_c_p, in_c_p);
00472       art::dotMultiplyDataData<double>(outDM_c_p, data_c_p, in_c_p);
00473       rst::subtract(out_c_p, outSM_c_p, outDM_c_p);
00474       if (rst::vectorNorm(&out_c_p[0], out_c_p.size(), NORM_ONE) > zero) {
00475         *outStream << "\n\nINCORRECT dotMultiplyDataData (1): check dot multiply for scalars vs. scalar multiply\n\n";
00476         errorFlag = -1000;
00477       }
00478       art::dotMultiplyDataData<double>(out_c_p, data_c_p_d, in_c_p_d);
00479       if (rst::vectorNorm(&out_c_p[0], out_c_p.size(), NORM_ONE) > zero) {
00480         *outStream << "\n\nINCORRECT dotMultiplyDataData (2): check dot multiply of orthogonal vectors\n\n";
00481         errorFlag = -1000;
00482       }
00483       art::dotMultiplyDataData<double>(out_c_p, data_c_p_d_d, in_c_p_d_d);
00484       if ((rst::vectorNorm(&out_c_p[0], out_c_p.size(), NORM_INF) - d1*d2) > zero) {
00485         *outStream << "\n\nINCORRECT dotMultiplyDataData (3): check dot multiply for tensors of 1s\n\n";
00486         errorFlag = -1000;
00487       }
00488       art::scalarMultiplyDataData<double>(outSM_c_p, data_c_1, in_c_p);
00489       art::dotMultiplyDataData<double>(outDM_c_p, data_c_1, in_c_p);
00490       rst::subtract(out_c_p, outSM_c_p, outDM_c_p);
00491       if (rst::vectorNorm(&out_c_p[0], out_c_p.size(), NORM_ONE) > zero) {
00492         *outStream << "\n\nINCORRECT dotMultiplyDataData (4): check dot multiply for scalars vs. scalar multiply\n\n";
00493         errorFlag = -1000;
00494       }
00495       art::dotMultiplyDataData<double>(out_c_p, data_c_1_d, in_c_p_d);
00496       if (rst::vectorNorm(&out_c_p[0], out_c_p.size(), NORM_ONE) > zero) {
00497         *outStream << "\n\nINCORRECT dotMultiplyDataData (5): check dot multiply of orthogonal vectors\n\n";
00498         errorFlag = -1000;
00499       }
00500       art::dotMultiplyDataData<double>(out_c_p, data_c_1_d_d, in_c_p_d_d);
00501       if ((rst::vectorNorm(&out_c_p[0], out_c_p.size(), NORM_INF) - d1*d2) > zero) {
00502         *outStream << "\n\nINCORRECT dotMultiplyDataData (6): check dot multiply for tensors of 1s\n\n";
00503         errorFlag = -1000;
00504       }
00505       } // end scope
00506 
00507 
00508       { // start scope
00509       *outStream << "\n************ Checking dotMultiplyDataData, (branch 2) ************\n";
00510 
00511       int c=5, p=9, d1=6, d2=14;
00512 
00513       FieldContainer<double> in_p(p);
00514       FieldContainer<double> in_p_d(p, d1);
00515       FieldContainer<double> in_p_d_d(p, d1, d2);
00516       FieldContainer<double> data_c_p(c, p);
00517       FieldContainer<double> data_c_p_d(c, p, d1);
00518       FieldContainer<double> data_c_p_d_d(c, p, d1, d2);
00519       FieldContainer<double> data_c_1(c, 1);
00520       FieldContainer<double> data_c_1_d(c, 1, d1);
00521       FieldContainer<double> data_c_1_d_d(c, 1, d1, d2);
00522       FieldContainer<double> out_c_p(c, p);
00523       FieldContainer<double> outSM_c_p(c, p);
00524       FieldContainer<double> outDM_c_p(c, p);
00525       double zero = INTREPID_TOL*10000.0;
00526 
00527       // fill with random numbers
00528       for (int i=0; i<in_p.size(); i++) {
00529         in_p[i] = Teuchos::ScalarTraits<double>::random();
00530       }
00531       // fill with alternating 1's and -1's
00532       for (int i=0; i<in_p_d.size(); i++) {
00533         in_p_d[i] = std::pow((double)-1.0, (int)i);
00534       }
00535       // fill with 1's
00536       for (int i=0; i<in_p_d_d.size(); i++) {
00537         in_p_d_d[i] = 1.0;
00538       }
00539       // fill with random numbers
00540       for (int i=0; i<data_c_p.size(); i++) {
00541         data_c_p[i] = Teuchos::ScalarTraits<double>::random();
00542       }
00543       // fill with 1's
00544       for (int i=0; i<data_c_1.size(); i++) {
00545         data_c_1[i] = Teuchos::ScalarTraits<double>::random();
00546       }
00547       // fill with 1's
00548       for (int i=0; i<data_c_p_d.size(); i++) {
00549         data_c_p_d[i] = 1.0;
00550       }
00551       // fill with 1's
00552       for (int i=0; i<data_c_1_d.size(); i++) {
00553         data_c_1_d[i] = 1.0;
00554       }
00555       // fill with 1's
00556       for (int i=0; i<data_c_p_d_d.size(); i++) {
00557         data_c_p_d_d[i] = 1.0;
00558       }
00559       // fill with 1's
00560       for (int i=0; i<data_c_1_d_d.size(); i++) {
00561         data_c_1_d_d[i] = 1.0;
00562       }
00563 
00564       art::scalarMultiplyDataData<double>(outSM_c_p, data_c_p, in_p);
00565       art::dotMultiplyDataData<double>(outDM_c_p, data_c_p, in_p);
00566       rst::subtract(out_c_p, outSM_c_p, outDM_c_p);
00567       if (rst::vectorNorm(&out_c_p[0], out_c_p.size(), NORM_ONE) > zero) {
00568         *outStream << "\n\nINCORRECT dotMultiplyDataData (7): check dot multiply for scalars vs. scalar multiply\n\n";
00569         errorFlag = -1000;
00570       }
00571       art::dotMultiplyDataData<double>(out_c_p, data_c_p_d, in_p_d);
00572       if (rst::vectorNorm(&out_c_p[0], out_c_p.size(), NORM_ONE) > zero) {
00573         *outStream << "\n\nINCORRECT dotMultiplyDataData (8): check dot multiply of orthogonal vectors\n\n";
00574         errorFlag = -1000;
00575       }
00576       art::dotMultiplyDataData<double>(out_c_p, data_c_p_d_d, in_p_d_d);
00577       if ((rst::vectorNorm(&out_c_p[0], out_c_p.size(), NORM_INF) - d1*d2) > zero) {
00578         *outStream << "\n\nINCORRECT dotMultiplyDataData (9): check dot multiply for tensors of 1s\n\n";
00579         errorFlag = -1000;
00580       }
00581       art::scalarMultiplyDataData<double>(outSM_c_p, data_c_1, in_p);
00582       art::dotMultiplyDataData<double>(outDM_c_p, data_c_1, in_p);
00583       rst::subtract(out_c_p, outSM_c_p, outDM_c_p);
00584       if (rst::vectorNorm(&out_c_p[0], out_c_p.size(), NORM_ONE) > zero) {
00585         *outStream << "\n\nINCORRECT dotMultiplyDataData (10): check dot multiply for scalars vs. scalar multiply\n\n";
00586         errorFlag = -1000;
00587       }
00588       art::dotMultiplyDataData<double>(out_c_p, data_c_1_d, in_p_d);
00589       if (rst::vectorNorm(&out_c_p[0], out_c_p.size(), NORM_ONE) > zero) {
00590         *outStream << "\n\nINCORRECT dotMultiplyDataData (11): check dot multiply of orthogonal vectors\n\n";
00591         errorFlag = -1000;
00592       }
00593       art::dotMultiplyDataData<double>(out_c_p, data_c_1_d_d, in_p_d_d);
00594       if ((rst::vectorNorm(&out_c_p[0], out_c_p.size(), NORM_INF) - d1*d2) > zero) {
00595         *outStream << "\n\nINCORRECT dotMultiplyDataData (12): check dot multiply for tensors of 1s\n\n";
00596         errorFlag = -1000;
00597       }
00598       } // end scope
00599 
00600       /******************************************/
00601       *outStream << "\n";
00602   }
00603   catch (std::logic_error err) {
00604     *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
00605     *outStream << err.what() << '\n';
00606     *outStream << "-------------------------------------------------------------------------------" << "\n\n";
00607     errorFlag = -1000;
00608   };
00609 
00610 
00611   if (errorFlag != 0)
00612     std::cout << "End Result: TEST FAILED\n";
00613   else
00614     std::cout << "End Result: TEST PASSED\n";
00615 
00616   // reset format state of std::cout
00617   std::cout.copyfmt(oldFormatState);
00618 
00619   return errorFlag;
00620 }