Thyra_VectorOpTester.hpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //    Thyra: Interfaces and Support for Abstract Numerical Algorithms
00005 //                 Copyright (2004) 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 Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef THYRA_VECTOROPTESTER_HPP
00030 #define THYRA_VECTOROPTESTER_HPP
00031 
00032 #include "Thyra_VectorImpl.hpp"
00033 #include "Thyra_VectorSpaceImpl.hpp"
00034 #include "Thyra_TesterBase.hpp"
00035 #include "Teuchos_ScalarTraits.hpp"
00036 #include "Teuchos_CommHelpers.hpp"
00037 #include "Thyra_ProductVectorSpaceBase.hpp"
00038 #include "Thyra_LinearCombinationImpl.hpp"
00039 
00040 namespace Thyra
00041 {
00042 using namespace Teuchos;
00043 using std::ostream;
00044 
00049 template <class Scalar>
00050 class VectorOpTester : public TesterBase<Scalar>
00051 {
00052 public:
00054   typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
00055 
00057   VectorOpTester(const RCP<const Comm<int> >& comm,
00058     const VectorSpace<Scalar>& space,
00059     Teuchos::RCP<Teuchos::FancyOStream>& out,
00060     const TestSpecifier<Scalar>& spec);
00061 
00063   bool runAllTests() const ;
00064 
00066   bool sumTest() const ;
00067 
00069   bool setElementUsingBracketTest() const ;
00070 
00072   bool dotStarTest() const ;
00073 
00075   bool scalarMultTest() const ;
00076 
00078   bool overloadedUpdateTest() const ;
00079 
00081   bool reciprocalTest() const ;
00082 
00084   bool minQuotientTest() const ;
00085 
00087   bool addScalarTest() const ;
00088 
00090   bool compareToScalarTest() const ;
00091 
00093   bool constraintMaskTest() const ;
00094 
00095 private:
00096   TestSpecifier<Scalar> spec_;
00097 };
00098 
00099 template <class Scalar> 
00100 inline VectorOpTester<Scalar>
00101 ::VectorOpTester(const RCP<const Comm<int> >& comm,
00102   const VectorSpace<Scalar>& space,
00103   Teuchos::RCP<Teuchos::FancyOStream>& out,
00104   const TestSpecifier<Scalar>& spec)
00105   : TesterBase<Scalar>(comm, space, 1, out),
00106     spec_(spec)
00107 {;}
00108 
00109 template <class Scalar> 
00110 inline bool VectorOpTester<Scalar>
00111 ::runAllTests() const
00112 {
00113   bool pass = true;
00114 
00115   pass = setElementUsingBracketTest() && pass;
00116   pass = sumTest() && pass;
00117   pass = dotStarTest() && pass;
00118   pass = scalarMultTest() && pass;
00119   pass = overloadedUpdateTest() && pass;
00120   pass = reciprocalTest() && pass;
00121   pass = minQuotientTest() && pass;
00122   pass = constraintMaskTest() && pass;
00123   pass = compareToScalarTest() && pass;
00124 
00125   return pass;
00126 }
00127 
00128 template <class Scalar> 
00129 inline bool VectorOpTester<Scalar>
00130 ::sumTest() const 
00131 {
00132 
00133   using std::endl;
00134 
00135   typedef Teuchos::ScalarTraits<Scalar> ST;
00136   ScalarMag err = ST::magnitude(ST::zero());
00137     
00138   this->out() << "running vector addition test..." << endl;
00139 
00140   Vector<Scalar> a = this->space().createMember();
00141   Vector<Scalar> b = this->space().createMember();
00142   Vector<Scalar> x = this->space().createMember();
00143   Vector<Scalar> y = this->space().createMember();
00144   randomizeVec(a);
00145   randomizeVec(b);
00146     
00147   /* do the operation with member functions */
00148   x = a + b ;
00149     
00150   /* do the operation elementwise. For off-processor elements, 
00151    * this is a no-op */
00152   for (int i=0; i<this->space().dim(); i++)
00153   {
00154     Scalar a_i = a[i];
00155     Scalar b_i = b[i];
00156     y[i] = a_i + b_i;
00157   }
00158   err = normInf(x-y);
00159   return this->checkTest(spec_, err, "vector addition");
00160 }
00161 
00162  
00163 
00164 
00165 template <class Scalar> 
00166 inline bool VectorOpTester<Scalar>
00167 ::setElementUsingBracketTest() const 
00168 {
00169 
00170   using std::endl;
00171 
00172   typedef Teuchos::ScalarTraits<Index> IST;
00173   typedef Teuchos::ScalarTraits<Scalar> ST;
00174 
00175   /* this test requires a comparison operator to make sense */
00176   if (ST::isComparable && spec_.doTest())
00177   {
00178     ScalarMag err = ST::magnitude(ST::zero());
00179 
00180     Vector<Scalar> a = this->space().createMember();
00181 
00182     /* We will load a vector with zeros, then set a randomly 
00183      * selected element to 1.0 and 
00184      * another to -1.0. Run minloc and maxloc to test that the minimum and
00185      * maximum values are at the expected locations. */
00186 
00187     zeroOut(a);
00188 
00189     int N = this->space().dim();
00190     int nTrials = 50;
00191 
00192     /* pick a seed in a way that is deterministic across processors */
00193     unsigned int seed = 107*N % 101;
00194     ST::seedrandom( seed );
00195 
00196     for (int t=0; t<nTrials; t++)
00197     {
00198       zeroOut(a);
00199       Index minIndex = IST::random() % N;
00200       Index maxIndex = IST::random() % N;
00201       /* skip cases where we've generated identical indices */
00202       if (minIndex == maxIndex) continue;
00203       a[minIndex] = -ST::one();
00204       a[maxIndex] = ST::one();
00205       Index findMin = -1;
00206       Index findMax = -1;
00207       Scalar minVal = Thyra::minloc(a, findMin);
00208       Scalar maxVal = Thyra::maxloc(a, findMax);
00209       err += ST::magnitude(-ST::one() - minVal);
00210       err += ST::magnitude(ST::one() - maxVal);
00211       err += ST::magnitude(findMin - minIndex);
00212       err += ST::magnitude(findMax - maxIndex);
00213     }
00214 
00215     return this->checkTest(spec_, err, "bracket operator");
00216   }
00217   else
00218   {
00219     this->out() << "skipping bracket operator test..." << endl;
00220     return true;
00221   }
00222 }
00223 
00224 /* specialize the set element test for complex types to a no-op. 
00225  * This is because minloc and maxloc do not exist for complex vectors 
00226  */
00227 #if defined(HAVE_COMPLEX) && defined(HAVE_TEUCHOS_COMPLEX)
00228 template <> 
00229 inline bool VectorOpTester<std::complex<double> >
00230 ::setElementUsingBracketTest() const 
00231 {
00232   using std::endl;
00233   this->out() << "skipping vector bracket operator test..." << endl;
00234   return true;
00235 }
00236 
00237 template <> 
00238 inline bool VectorOpTester<std::complex<float> >
00239 ::setElementUsingBracketTest() const 
00240 {
00241   using std::endl;
00242   this->out() << "skipping vector bracket operator test..." << endl;
00243   return true;
00244 }
00245 #endif
00246 
00247 
00248   
00249 
00250 template <class Scalar> 
00251 inline bool VectorOpTester<Scalar>
00252 ::dotStarTest() const 
00253 {
00254 
00255   using std::endl;
00256 
00257   typedef Teuchos::ScalarTraits<Scalar> ST;
00258   ScalarMag err = ST::magnitude(ST::zero());
00259 
00260   this->out() << "running vector dotStar test..." << endl;
00261     
00262   Vector<Scalar> a = this->space().createMember();
00263   Vector<Scalar> b = this->space().createMember();
00264   Vector<Scalar> x = this->space().createMember();
00265   Vector<Scalar> y = this->space().createMember();
00266   Vector<Scalar> z = this->space().createMember();
00267   randomizeVec(a);
00268   randomizeVec(b);
00269     
00270   x = dotStar(a,b);
00271   z = dotSlash(x, b);
00272     
00273   err = normInf(a-z);
00274     
00275   this->out() << "|dotStar error|=" << err << endl;
00276   return this->checkTest(spec_, err, "dot-star");
00277 }
00278 
00279 
00280 
00281 template <class Scalar> 
00282 inline bool VectorOpTester<Scalar>
00283 ::scalarMultTest() const 
00284 {
00285 
00286   using std::endl;
00287 
00288   typedef Teuchos::ScalarTraits<Scalar> ST;
00289   ScalarMag err = ST::magnitude(ST::zero());
00290 
00291   this->out() << "running vector scalarMult test..." << endl;
00292     
00293   Vector<Scalar> a = this->space().createMember();
00294   Vector<Scalar> x = this->space().createMember();
00295   Vector<Scalar> y = this->space().createMember();
00296   randomizeVec(a);
00297     
00298     
00299   /* do the operation with member functions */
00300   Scalar alpha = 3.14;
00301   x = alpha * a;
00302     
00303   /* do the operation elementwise */
00304     
00305   for (int i=0; i<this->space().dim(); i++)
00306   {
00307     Scalar a_i = a[i];
00308     y[i]= alpha *a_i;
00309   }
00310     
00311   err = normInf(x-y);
00312     
00313   this->out() << "|scalarMult error|=" << err << endl;
00314   return this->checkTest(spec_, err, "scalar multiplication");
00315 
00316 }
00317  
00318 template <class Scalar> 
00319 inline bool VectorOpTester<Scalar>
00320 ::overloadedUpdateTest() const 
00321 {
00322 
00323   using std::endl;
00324 
00325   typedef Teuchos::ScalarTraits<Scalar> ST;
00326   ScalarMag err = ST::magnitude(ST::zero());
00327 
00328   this->out() << "running vector overloadedUpdate test..." << endl;
00329     
00330   Vector<Scalar> a = this->space().createMember();
00331   Vector<Scalar> b = this->space().createMember();
00332   Vector<Scalar> x = this->space().createMember();
00333   Vector<Scalar> y = this->space().createMember();
00334   randomizeVec(a);
00335   randomizeVec(b);
00336     
00337   /* do the operation with member functions */
00338   Scalar alpha = 3.14;
00339   Scalar beta = 1.4;
00340   x = alpha*a + beta*b;
00341     
00342   /* do the operation elementwise */
00343   for (int i=0; i<this->space().dim(); i++)
00344   {
00345     Scalar a_i = a[i];
00346     Scalar b_i = b[i];
00347     y[i] = alpha*a_i + beta*b_i;
00348   }
00349     
00350   err = normInf(x-y);
00351     
00352   this->out() << "|overloadedUpdate error|=" << err << endl;
00353   return this->checkTest(spec_, err, "overloaded update");
00354 
00355 }
00356 
00357 template <class Scalar> 
00358 inline bool VectorOpTester<Scalar>
00359 ::reciprocalTest() const 
00360 {
00361 
00362   using std::endl;
00363 
00364   typedef Teuchos::ScalarTraits<Scalar> ST;
00365   ScalarMag err = ST::magnitude(ST::zero());
00366 
00367   this->out() << "running vector reciprocal test..." << endl;
00368 
00369   Vector<Scalar> a = this->space().createMember();
00370   randomizeVec(a);
00371     
00372   Vector<Scalar> y = this->space().createMember();
00373 
00374   /* load the operation elementwise */
00375   int localDenomsAreOK = true;
00376   for (int i=0; i<this->space().dim(); i++)
00377   {
00378     if (!indexIsLocal(this->space(), i)) continue;
00379     Scalar a_i = a[i];
00380     if (a_i != Teuchos::ScalarTraits<Scalar>::zero()) 
00381     {
00382       y[i] = Teuchos::ScalarTraits<Scalar>::one()/a_i;
00383     }
00384     else
00385     {
00386       localDenomsAreOK=false;
00387       y[i] = a_i;
00388     }
00389   }
00390         
00391   Vector<Scalar> x = this->space().createMember();
00392   int tDenomsAreOK = Thyra::VInvTest(*(a.constPtr()), x.ptr().get());
00393   err = norm2(x - y);
00394 
00395   int allDenomsAreOK;
00396     
00397   reduceAll(this->comm(), Teuchos::REDUCE_AND, 1, 
00398     &localDenomsAreOK, &allDenomsAreOK);
00399 
00400   this->out() << "|reciprocal error|=" << err << endl;
00401   this->out() << "denom check test=" << (allDenomsAreOK != tDenomsAreOK)
00402               << endl;
00403   return this->checkTest(spec_, err + (allDenomsAreOK != tDenomsAreOK), 
00404     "reciprocal");
00405 
00406 }
00407 
00408 template <class Scalar> 
00409 inline bool VectorOpTester<Scalar>
00410 ::minQuotientTest() const 
00411 {
00412 
00413   using std::endl;
00414 
00415   typedef Teuchos::ScalarTraits<Scalar> ST;
00416   ScalarMag err = ST::magnitude(ST::zero());
00417 
00418   this->out() << "running vector minQuotient test..." << endl;
00419     
00420   Vector<Scalar> a = this->space().createMember();
00421   Vector<Scalar> b = this->space().createMember();
00422     
00423   randomizeVec(a);
00424   randomizeVec(b);
00425     
00426   /* perform the operation elementwise */
00427 
00428   Scalar minQLocal = Teuchos::ScalarTraits<Scalar>::rmax();
00429   for (int i=0; i<this->space().dim(); i++)
00430   {
00431     if (!indexIsLocal(this->space(), i)) continue;
00432     Scalar a_i = a[i];
00433     Scalar b_i = b[i];
00434     if (b_i != Teuchos::ScalarTraits<Scalar>::zero())
00435     {
00436       Scalar q = a_i/b_i;
00437       if (q < minQLocal) minQLocal = q;
00438     }
00439   }
00440 
00441   Scalar minQ = minQLocal;
00442     
00443   reduceAll(this->comm(), Teuchos::REDUCE_MIN, 1, 
00444     &minQLocal, &minQ);
00445 
00446   Scalar tMinQ = Thyra::VMinQuotient(*(a.constPtr()), *(b.constPtr()));
00447 
00448   this->out() << "trilinos minQ = " << tMinQ << endl;
00449   this->out() << "elemwise minQ = " << minQ << endl;
00450 
00451   err = ST::magnitude(minQ - tMinQ);
00452     
00453   this->out() << "min quotient error=" << err << endl;
00454 
00455   return this->checkTest(spec_, err, "min quotient");
00456 
00457 }
00458 
00459 /* specialize the min quotient test for complex types to a no-op. 
00460  * This is because min function does not exist for complex vectors 
00461  */
00462 #if defined(HAVE_COMPLEX) && defined(HAVE_TEUCHOS_COMPLEX)
00463 template <> 
00464 inline bool VectorOpTester<std::complex<double> >
00465 ::minQuotientTest() const 
00466 {
00467   using std::endl;
00468   this->out() << "skipping min quotient test..." << endl;
00469   return true;
00470 }
00471 
00472 template <> 
00473 inline bool VectorOpTester<std::complex<float> >
00474 ::minQuotientTest() const 
00475 {
00476   using std::endl;
00477   this->out() << "skipping min quotient test..." << endl;
00478   return true;
00479 }
00480 #endif
00481 
00482 template <class Scalar> 
00483 inline bool VectorOpTester<Scalar>
00484 ::constraintMaskTest() const 
00485 {
00486 
00487   using std::endl;
00488 
00489   typedef Teuchos::ScalarTraits<Scalar> ST;
00490   ScalarMag err = ST::magnitude(ST::zero());
00491 
00492   this->out() << "running vector constraintMask test..." << endl;
00493 
00494   Vector<Scalar> a = this->space().createMember();
00495   Vector<Scalar> c = this->space().createMember();
00496   randomizeVec(a);
00497     
00498   Vector<Scalar> y = this->space().createMember();
00499   Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
00500     
00501   /* do the operation elementwise */
00502   int allFeasible = true;
00503   for (int i=0; i<this->space().dim(); i++)
00504   {
00505     if (!indexIsLocal(this->space(), i)) continue;
00506     int feasible = true;
00507     Scalar a_i = a[i];
00508     switch(i%4)
00509     {
00510       case 0:
00511         c[i] = -2.0;
00512         feasible = a_i < zero;
00513         break;
00514       case 1:
00515         c[i] = -1.0;
00516         feasible = a_i <= zero;
00517         break;
00518       case 2:
00519         c[i] = 1.0;
00520         feasible = a_i > zero;
00521         break;
00522       case 3:
00523         c[i] = 2.0;
00524         feasible = a_i >= zero;
00525         break;
00526       default:
00527         TEST_FOR_EXCEPTION(true, std::logic_error, "impossible!");
00528     }
00529     y[i] = (Scalar) !feasible;
00530     allFeasible = allFeasible && feasible;
00531   }
00532     
00533   Vector<Scalar> m = this->space().createMember();
00534   int tAllFeasible = Thyra::VConstrMask(*(a.constPtr()), *(c.constPtr()), m.ptr().get());
00535   err = norm2(m - y);
00536     
00537   int localAllFeas = allFeasible;
00538   this->out() << "local all feas=" << localAllFeas << endl;
00539     
00540   reduceAll(this->comm(), Teuchos::REDUCE_AND, 1, 
00541     &localAllFeas, &allFeasible);
00542   this->out() << "global all feas=" << allFeasible << endl;
00543     
00544   this->out() << "|constraintMask error|=" << err << endl;
00545   if (err > spec_.errorTol())
00546   {
00547     this->out() << "vector constraintMask test FAILED: tol = " 
00548                 << spec_.errorTol() << endl;
00549     return false;
00550   }
00551   else if (allFeasible != tAllFeasible)
00552   {
00553     this->out() << "vector constraintMask test FAILED: trilinosFeas="
00554                 << tAllFeasible << ", feas=" << allFeasible << endl;
00555     return false;
00556   }
00557   else if (err > spec_.warningTol())
00558   {
00559     this->out() << "WARNING: vector constraintMask test could not beat tol = " 
00560                 << spec_.warningTol() << endl;
00561   }
00562   else
00563   {
00564     this->out() << "vector constraintMask test PASSED: error=" << err << ", tol = " 
00565                 << spec_.errorTol() << endl;
00566   }
00567   return true;
00568 }
00569   
00570 
00571 /* specialize the constraint mask  test for complex types to a no-op. 
00572  * This is because min function does not exist for complex vectors 
00573  */
00574 #if defined(HAVE_COMPLEX) && defined(HAVE_TEUCHOS_COMPLEX)
00575 template <> 
00576 inline bool VectorOpTester<std::complex<double> >
00577 ::constraintMaskTest() const 
00578 {
00579   using std::endl;
00580   this->out() << "skipping constraint mask test..." << endl;
00581   return true;
00582 }
00583 
00584 template <> 
00585 inline bool VectorOpTester<std::complex<float> >
00586 ::constraintMaskTest() const 
00587 {
00588   using std::endl;
00589   this->out() << "skipping constraint mask test..." << endl;
00590   return true;
00591 }
00592 
00593 #endif
00594 
00595 template <class Scalar> 
00596 inline bool VectorOpTester<Scalar>
00597 ::compareToScalarTest() const 
00598 {
00599 
00600   using std::endl;
00601 
00602   typedef Teuchos::ScalarTraits<Scalar> ST;
00603   ScalarMag err = ST::magnitude(ST::zero());
00604 
00605   this->out() << "running vector compare-to-scalar test..." << endl;
00606     
00607   Vector<Scalar> a = this->space().createMember();
00608   Vector<Scalar> x = this->space().createMember();
00609   Vector<Scalar> y = this->space().createMember();
00610   randomizeVec(a);
00611     
00612   /* do the operation with member functions */
00613   Scalar s = 0.5;
00614   Thyra::VCompare(s, *(a.constPtr()), x.ptr().get());
00615     
00616   /* do the operation elementwise */    
00617   for (int i=0; i<this->space().dim(); i++)
00618   {
00619     if (!indexIsLocal(this->space(), i)) continue;
00620     Scalar a_i = a[i];
00621     y[i] = std::fabs(a_i) >= s ;
00622   }
00623     
00624   err = normInf(x-y);
00625 
00626   this->out() << "|compare-to-scalar error|=" << err << endl;
00627   return this->checkTest(spec_, err, "compare-to-scalar");
00628 
00629 }
00630 
00631 
00632 /* specialize the compare to scalar test for complex types to a no-op. 
00633  * This is because comparison function does not exist for complex vectors 
00634  */
00635 #if defined(HAVE_COMPLEX) && defined(HAVE_TEUCHOS_COMPLEX)
00636 template <> 
00637 inline bool VectorOpTester<std::complex<double> >
00638 ::compareToScalarTest() const 
00639 {
00640   using std::endl;
00641   this->out() << "skipping compare-to-scalar test..." << endl;
00642   return true;
00643 }
00644 
00645 template <> 
00646 inline bool VectorOpTester<std::complex<float> >
00647 ::compareToScalarTest() const 
00648 {
00649   using std::endl;
00650   this->out() << "skipping compare-to-scalar test..." << endl;
00651   return true;
00652 }
00653 #endif
00654 
00655 
00656 }
00657 
00658 #endif

Generated on Tue Oct 20 12:47:27 2009 for Thyra Package Browser (Single Doxygen Collection) by doxygen 1.4.7