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

Generated on Thu Sep 18 12:33:03 2008 for Thyra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1