Thyra Package Browser (Single Doxygen Collection): Thyra_LinearCombinationTester.hpp Source 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_LINEARCOMBINATIONTESTER_HPP
00030 #define THYRA_LINEARCOMBINATIONTESTER_HPP
00031 
00032 #include "Thyra_LinearOperatorImpl.hpp"
00033 #include "Thyra_LinearCombinationImpl.hpp"
00034 #include "Thyra_VectorImpl.hpp"
00035 #include "Thyra_TestSpecifier.hpp"
00036 #include "Thyra_TesterBase.hpp"
00037 #include "Teuchos_ScalarTraits.hpp"
00038 
00039 #define TESTER(form1, form2)\
00040   {\
00041     Vector<Scalar> _val1 = form1;\
00042     TEST_FOR_EXCEPT(_val1.constPtr().get()==0);\
00043     Vector<Scalar> _val2 = form2;\
00044     TEST_FOR_EXCEPT(_val2.constPtr().get()==0);\
00045     ScalarMag err = norm2(_val1-_val2);\
00046     if (!checkTest(spec_, err, "[" #form1 "] == [" #form2 "]")) pass = false;\
00047   }
00048 
00049 namespace Thyra
00050 {
00051   using Teuchos::RefCountPtr;
00052   using Teuchos::ScalarTraits;
00053 
00055   template <class Scalar>
00056   class LinearCombinationTester : public TesterBase<Scalar>
00057   {
00058   public:
00059 
00061     LinearCombinationTester(const RefCountPtr<const Comm<int> >& comm,
00062                             const VectorSpace<Scalar>& vecSpace, 
00063                             Teuchos::RefCountPtr<Teuchos::FancyOStream>& out,
00064                             const TestSpecifier<Scalar>& spec);
00065 
00067     bool runAllTests() const ;
00068 
00069 
00070   private:
00071 
00073     bool nonModifyingOpTests() const ;
00074 
00076     bool selfModifyingOpTests() const ;
00077 
00078     TestSpecifier<Scalar> spec_;
00079 
00080   };
00081 
00082   template <class Scalar> 
00083   inline LinearCombinationTester<Scalar>
00084   ::LinearCombinationTester(const RefCountPtr<const Comm<int> >& comm,
00085                             const VectorSpace<Scalar>& vecSpace,
00086                             Teuchos::RefCountPtr<Teuchos::FancyOStream>& out,
00087                             const TestSpecifier<Scalar>& spec)
00088     : TesterBase<Scalar>(comm, vecSpace, vecSpace.dim(), out), spec_(spec)
00089   {;}
00090 
00091   template <class Scalar> 
00092   inline bool LinearCombinationTester<Scalar>
00093   ::runAllTests() const
00094   {
00095     bool pass = true;
00096 
00097     pass = this->nonModifyingOpTests() && pass;
00098     pass = this->selfModifyingOpTests() && pass;
00099 
00100     return pass;
00101   }
00102 
00103   template <class Scalar> 
00104   inline bool LinearCombinationTester<Scalar>
00105   ::nonModifyingOpTests() const
00106   {
00107     bool pass = true;
00108     typedef typename Teuchos::ScalarTraits<Scalar> ST;
00109     typedef typename ST::magnitudeType ScalarMag;
00110 
00111     LinearOperator<Scalar> A = this->randomDenseOp();
00112     LinearOperator<Scalar> B = this->randomDenseOp();
00113     LinearOperator<Scalar> C = this->randomDenseOp();
00114 
00115     Vector<Scalar> x = A.domain().createMember();
00116     Vector<Scalar> y = A.domain().createMember();
00117     Vector<Scalar> z = A.domain().createMember();
00118 
00119     randomizeVec(x);
00120     randomizeVec(y);
00121     randomizeVec(z);
00122 
00123     this->out() << "starting linear combination tests" << endl;
00124 
00125     TESTER(x*Scalar(2.0), Scalar(2.0)*x);
00126 
00127     TESTER(Scalar(2.0)*(x + y), Scalar(2.0)*x + Scalar(2.0)*y);
00128 
00129     TESTER(Scalar(2.0)*(x - y), Scalar(2.0)*x - Scalar(2.0)*y);
00130 
00131     TESTER((x + y)*Scalar(2.0), Scalar(2.0)*x + Scalar(2.0)*y);
00132 
00133     TESTER((x - y)*Scalar(2.0), Scalar(2.0)*x - Scalar(2.0)*y);
00134 
00135     TESTER(Scalar(2.0)*(x - y), -Scalar(2.0)*(y - x));
00136 
00137     TESTER(Scalar(0.25)*(Scalar(2.0)*(x + y) - Scalar(2.0)*(x - y)), y);
00138 
00139     TESTER((Scalar(2.0)*A)*x, Scalar(2.0)*(A*x));
00140 
00141     TESTER(Scalar(2.0)*(A*x), (A*x)*Scalar(2.0));
00142 
00143     TESTER(A*(B*x), (A*B)*x);
00144 
00145     TESTER(Scalar(2.0)*A*(B*x), A*(B*(Scalar(2.0)*x)));
00146 
00147     TESTER(Scalar(3.0)*(Scalar(2.0)*A)*x, Scalar(6.0)*(A*x));
00148 
00149     TESTER(A*x + y, y + A*x);
00150 
00151     TESTER(z + (A*x + B*y), (B*y + A*x) + z);
00152 
00153     TESTER(z - (A*x + B*y), Scalar(-1.0)*((B*y + A*x) - z));
00154 
00155     TESTER(C*z + (A*x + B*y), (B*y + A*x) + C*z);
00156 
00157     TESTER(C*z - (A*x + B*y), Scalar(-1.0)*((B*y + A*x) - C*z));
00158 
00159     TESTER(Scalar(2.0)*z + (A*x + B*y), (B*y + A*x) + Scalar(2.0)*z);
00160 
00161     TESTER(Scalar(2.0)*z - (A*x + B*y), Scalar(-1.0)*((B*y + A*x) - Scalar(2.0)*z));
00162 
00163     TESTER(A*x - y, Scalar(-1.0)*(y - A*x));
00164 
00165     TESTER(A*x + B*y, B*y + A*x);
00166 
00167     TESTER(A*x - B*y - A*x + B*y +z, z);
00168 
00169     TESTER(Scalar(2.0)*(A*x + y), Scalar(2.0)*A*x + Scalar(2.0)*y);
00170 
00171     TESTER(Scalar(2.0)*(A*x + B*y), A*x + B*y + A*x + B*y);
00172 
00173     TESTER(Scalar(2.0)*(y + A*x), Scalar(2.0)*y + Scalar(2.0)*(A*x));
00174 
00175     TESTER(x + Scalar(2.0)*A*y, x + Scalar(2.0)*(A*y));
00176 
00177     TESTER(Scalar(2.0)*A*y + x, Scalar(2.0)*(A*y) + x);
00178 
00179     TESTER(Scalar(2.0)*(A*x + B*y), Scalar(2.0)*A*x + Scalar(2.0)*B*y);
00180 
00181     TESTER(Scalar(2.0)*(A*x - B*y), Scalar(2.0)*A*x - Scalar(2.0)*B*y);
00182 
00183     TESTER(Scalar(2.0)*(A*x + B*y + z), Scalar(2.0)*A*x + Scalar(2.0)*B*y + Scalar(2.0)*z);
00184 
00185     TESTER(Scalar(2.0)*(A*x + Scalar(3.0)*B*y), Scalar(2.0)*A*x + Scalar(6.0)*B*y);
00186 
00187     TESTER(Scalar(2.0)*(A*x + Scalar(3.0)*(z + B*y)), Scalar(2.0)*A*x + Scalar(6.0)*B*y + Scalar(6.0)*z);
00188 
00189     TESTER(Scalar(2.0)*(z + A*x + B*y + z), Scalar(2.0)*A*x + Scalar(2.0)*B*y + Scalar(4.0)*z);
00190 
00191     TESTER(Scalar(2.0)*(Scalar(3.0)*(z + A*x) + B*y), Scalar(6.0)*z + Scalar(6.0)*A*x + Scalar(2.0)*B*y);
00192 
00193     TESTER(Scalar(2.0)*(Scalar(3.0)*(z + A*x) + Scalar(4.0)*(B*y + z)), Scalar(6.0)*z + Scalar(6.0)*A*x + Scalar(8.0)*B*y + Scalar(8.0)*z);
00194     
00195     TESTER((A*x + B*y) + (A*y + B*x), (A + B)*x + (A+B)*y);
00196     TESTER((A*x + B*y) - (A*y + B*x), A*x - A*y + B*y - B*x);
00197 
00198     TESTER((A*x + B*y) + Scalar(2.0)*(A*y + B*x), A*(x + Scalar(2.0)*y) + B*(Scalar(2.0)*x + y));
00199     TESTER((A*x + B*y) - Scalar(2.0)*(A*y + B*x), A*(x - Scalar(2.0)*y) + B*(y - Scalar(2.0)*x));
00200 
00201     return pass;
00202   }
00203 
00204   template <class Scalar> 
00205   inline bool LinearCombinationTester<Scalar>
00206   ::selfModifyingOpTests() const
00207   {
00208     bool pass = true;
00209     typedef typename Teuchos::ScalarTraits<Scalar> ST;
00210     typedef typename ST::magnitudeType ScalarMag;
00211 
00212     LinearOperator<Scalar> A = this->randomDenseOp();
00213     LinearOperator<Scalar> B = this->randomDenseOp();
00214     LinearOperator<Scalar> C = this->randomDenseOp();
00215 
00216     Vector<Scalar> x = A.domain().createMember();
00217     Vector<Scalar> y = A.domain().createMember();
00218     Vector<Scalar> z = A.domain().createMember();
00219 
00220     randomizeVec(x);
00221     randomizeVec(y);
00222     randomizeVec(z);
00223 
00224     Vector<Scalar> a = copy(x);
00225     Vector<Scalar> b = copy(y);
00226     Vector<Scalar> c = copy(z);
00227     
00228     TestSpecifier<Scalar>
00229       specLooser(spec_.doTest(),1e+1*spec_.warningTol(),1e+2*spec_.errorTol());
00230 
00231     this->out() << "starting linear combination tests" << endl;
00232 
00233     x = Scalar(2.0)*A*x;
00234     ScalarMag err = norm2(x - Scalar(2.0)*A*a);
00235     if (!checkTest(spec_, err, "x=Scalar(2.0)*A*x")) pass = false;
00236 
00237     a = copy(x);
00238     x = x + y;
00239     err = norm2(x - (a + y));
00240     if (!checkTest(spec_, err, "x=x+y")) pass = false;
00241 
00242     a = copy(x);
00243     x = y + x;
00244     err = norm2(x - (y + a));
00245     if (!checkTest(spec_, err, "x=y+x")) pass = false;
00246 
00247     a = copy(x);
00248     x = A*x + B*y;
00249     err = norm2(x - (A*a + B*y));
00250     if (!checkTest(spec_, err, "x=A*x+B*y")) pass = false;
00251 
00252     a = copy(x);
00253     x = B*y + A*x;
00254     err = norm2(x - (B*y + A*a));
00255     if (!checkTest(spec_, err, "x=B*y+A*x")) pass = false;
00256 
00257     a = copy(x);
00258     x = A*x + (B*y + C*x);
00259     err = norm2(x - (A*a + (B*y + C*a)));
00260     if (!checkTest(specLooser, err, "x=A*x + (B*y + C*x)")) pass = false;
00261 
00262     a = copy(x);
00263     x = (A*x + B*y) + C*x;
00264     err = norm2(x - ((A*a + B*y) + C*a));
00265     if (!checkTest(specLooser, err, "x=(A*x + B*y) + C*x")) pass = false;
00266 
00267     /* test assignment of OpTimesLC into empty and non-empty vectors */
00268     Vector<Scalar> u;
00269     u = Scalar(2.0)*A*B*x;
00270     err = norm2(u - Scalar(2.0)*A*B*x);
00271 
00272     u = Scalar(2.0)*A*x;
00273     err = norm2(u - Scalar(2.0)*A*B*x);
00274 
00275     /* test assignment of LC2 into empty and non-empty vectors */
00276     Vector<Scalar> v;
00277     v = Scalar(2.0)*x + Scalar(3.0)*y;
00278     err = norm2(v - (Scalar(2.0)*x + Scalar(3.0)*y));
00279 
00280     v = Scalar(2.0)*x + Scalar(3.0)*y;
00281     err = norm2(v - (Scalar(2.0)*x + Scalar(3.0)*y));
00282 
00283     return pass;
00284   }
00285  
00286 }
00287 
00288 // We had better not leave a global macro named TESTSER just lying around
00289 // awating disaster!
00290 #undef TESTER
00291 
00292 #endif