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