00001
00002
00003
00004
00005
00006 #ifndef PLAYA_LINEARCOMBINATIONTESTER_HPP
00007 #define PLAYA_LINEARCOMBINATIONTESTER_HPP
00008
00009 #include "PlayaLinearOperatorDecl.hpp"
00010 #include "PlayaLinearCombinationImpl.hpp"
00011 #include "PlayaSimpleComposedOpDecl.hpp"
00012 #include "PlayaSimpleScaledOpDecl.hpp"
00013 #include "PlayaSimpleAddedOpDecl.hpp"
00014 #include "PlayaSimpleIdentityOpDecl.hpp"
00015 #include "PlayaOut.hpp"
00016 #include "PlayaVectorDecl.hpp"
00017 #include "PlayaTesterBase.hpp"
00018 #include "PlayaRandomSparseMatrixBuilderDecl.hpp"
00019 #include "PlayaTestSpecifier.hpp"
00020 #include "Teuchos_ScalarTraits.hpp"
00021
00022 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00023 #include "PlayaVectorImpl.hpp"
00024 #include "PlayaLinearOperatorImpl.hpp"
00025 #include "PlayaRandomSparseMatrixBuilderImpl.hpp"
00026 #include "PlayaSimpleComposedOpImpl.hpp"
00027 #include "PlayaSimpleScaledOpImpl.hpp"
00028 #include "PlayaSimpleAddedOpImpl.hpp"
00029 #include "PlayaSimpleTransposedOpImpl.hpp"
00030 #include "PlayaSimpleIdentityOpDecl.hpp"
00031 #endif
00032
00033 using namespace Playa;
00034 using namespace Teuchos;
00035
00036
00037
00038
00039
00040 #define TESTER(form1, form2)\
00041 {\
00042 Out::os() << "testing " #form1 << std::endl;\
00043 Vector<Scalar> _val1 = form1;\
00044 Out::os() << "testing " #form2 << std::endl;\
00045 Vector<Scalar> _val2 = form2;\
00046 Out::os() << "done testing... checking error" << std::endl;\
00047 ScalarMag err = (_val1-_val2).norm2();\
00048 if (!this->checkTest(spec_, err, "[" #form1 "] == [" #form2 "]")) pass = false;\
00049 }
00050
00051
00052
00053
00054
00055
00056
00057 namespace Playa
00058 {
00059
00060
00061 template <class Scalar>
00062 class LinearCombinationTester : public TesterBase<Scalar>
00063 {
00064 public:
00065
00066 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
00067
00068
00069 LinearCombinationTester(int nLocalRows, double onProcDensity,
00070 double offProcDensity,
00071 const VectorType<Scalar>& vecType,
00072 const TestSpecifier<Scalar>& spec);
00073
00074
00075 bool runAllTests() const ;
00076
00077
00078 private:
00079
00080
00081 bool nonModifyingOpTests() const ;
00082
00083
00084 bool selfModifyingOpTests() const ;
00085
00086 TestSpecifier<Scalar> spec_;
00087
00088 int nLocalRows_;
00089 double onProcDensity_;
00090 double offProcDensity_;
00091 VectorType<Scalar> vecType_;
00092
00093 };
00094
00095 template <class Scalar>
00096 inline LinearCombinationTester<Scalar>
00097 ::LinearCombinationTester(int nLocalRows, double onProcDensity,
00098 double offProcDensity,
00099 const VectorType<Scalar>& vecType,
00100 const TestSpecifier<Scalar>& spec)
00101 : TesterBase<Scalar>(),
00102 spec_(spec),
00103 nLocalRows_(nLocalRows),
00104 onProcDensity_(onProcDensity),
00105 offProcDensity_(offProcDensity),
00106 vecType_(vecType)
00107 {;}
00108
00109 template <class Scalar>
00110 inline bool LinearCombinationTester<Scalar>
00111 ::runAllTests() const
00112 {
00113 bool pass = true;
00114
00115 pass = nonModifyingOpTests() && pass;
00116 pass = selfModifyingOpTests() && pass;
00117
00118 return pass;
00119 }
00120
00121 template <class Scalar>
00122 inline bool LinearCombinationTester<Scalar>
00123 ::nonModifyingOpTests() const
00124 {
00125 bool pass = true;
00126
00127 RandomSparseMatrixBuilder<double> ABuilder(nLocalRows_, nLocalRows_,
00128 onProcDensity_, offProcDensity_,
00129 vecType_);
00130 RandomSparseMatrixBuilder<double> BBuilder(nLocalRows_, nLocalRows_,
00131 onProcDensity_, offProcDensity_,
00132 vecType_);
00133 RandomSparseMatrixBuilder<double> CBuilder(nLocalRows_, nLocalRows_,
00134 onProcDensity_, offProcDensity_,
00135 vecType_);
00136
00137
00138 LinearOperator<double> A = ABuilder.getOp();
00139 LinearOperator<Scalar> B = BBuilder.getOp();
00140 LinearOperator<Scalar> C = CBuilder.getOp();
00141
00142 Vector<Scalar> x = A.domain().createMember();
00143 Vector<Scalar> y = A.domain().createMember();
00144 Vector<Scalar> z = A.domain().createMember();
00145
00146 this->randomizeVec(x);
00147 this->randomizeVec(y);
00148 this->randomizeVec(z);
00149
00150
00151 TESTER(x*2.0, 2.0*x);
00152
00153 TESTER(2.0*(x + y), 2.0*x + 2.0*y);
00154
00155 TESTER(2.0*(x - y), 2.0*x - 2.0*y);
00156
00157 TESTER((2.0*x + y) - y, 2.0*x);
00158
00159 TESTER(-1.0*y + (2.0*x + y), 2.0*x);
00160
00161 TESTER((x + y)*2.0, 2.0*x + 2.0*y);
00162
00163 TESTER((x - y)*2.0, 2.0*x - 2.0*y);
00164
00165 TESTER(2.0*(x - y), -2.0*(y - x));
00166
00167 TESTER(0.25*(2.0*(x + y) - 2.0*(x - y)), y);
00168
00169 TESTER((2.0*A)*x, 2.0*(A*x));
00170
00171 TESTER(2.0*(A*x), (A*x)*2.0);
00172
00173 TESTER(A*(B*x), (A*B)*x);
00174
00175 TESTER(2.0*A*(B*x), A*(B*(2.0*x)));
00176
00177 TESTER(3.0*(2.0*A)*x, 6.0*(A*x));
00178
00179 TESTER(y + A*x, A*x + y);
00180
00181
00182 TESTER(z + (A*x + B*y), (B*y + A*x) + z);
00183
00184
00185 TESTER(z - (A*x + B*y), -1.0*((B*y + A*x) - z));
00186
00187 TESTER(C*z + (A*x + B*y), (B*y + A*x) + C*z);
00188
00189 TESTER(C*z - (A*x + B*y), -1.0*((B*y + A*x) - C*z));
00190
00191 TESTER(2.0*z + (A*x + B*y), (B*y + A*x) + 2.0*z);
00192
00193 TESTER(2.0*z - (A*x + B*y), -1.0*((B*y + A*x) - 2.0*z));
00194
00195 TESTER(A*x - y, -1.0*(y - A*x));
00196
00197 TESTER(A*x + B*y, B*y + A*x);
00198
00199 TESTER(A*x - B*y - A*x + B*y +z, z);
00200
00201 TESTER(2.0*(A*x + y), 2.0*A*x + 2.0*y);
00202
00203 TESTER(2.0*(A*x + B*y), A*x + B*y + A*x + B*y);
00204
00205 TESTER(2.0*(y + A*x), 2.0*y + 2.0*(A*x));
00206
00207 TESTER(x + 2.0*A*y, x + 2.0*(A*y));
00208
00209 TESTER(2.0*A*y + x, 2.0*(A*y) + x);
00210
00211 TESTER(2.0*A*(3.0*B)*y, 6.0*(A*B*y));
00212
00213 TESTER(2.0*A*(3.0*B)*y, 6.0*(A*(B*y)));
00214
00215 TESTER(2.0*A*(3.0*B + 2.0*A)*x, 6.0*A*B*x + 4.0*A*A*x );
00216
00217 TESTER(2.0*(A*x + B*y), 2.0*A*x + 2.0*B*y);
00218
00219 TESTER(2.0*(A*x - B*y), 2.0*A*x - 2.0*B*y);
00220
00221 TESTER(2.0*(A*x + B*y + z), 2.0*A*x + 2.0*B*y + 2.0*z);
00222
00223 TESTER(2.0*(A*x + 3.0*B*y), 2.0*A*x + 6.0*B*y);
00224
00225 TESTER(2.0*(A*x + 3.0*(z + B*y)), 2.0*A*x + 6.0*B*y + 6.0*z);
00226
00227 TESTER(2.0*(z + A*x + B*y + z), 2.0*A*x + 2.0*B*y + 4.0*z);
00228
00229 TESTER(2.0*(3.0*(z + A*x) + B*y), 6.0*z + 6.0*A*x + 2.0*B*y);
00230
00231 TESTER(2.0*(3.0*(z + A*x) + 4.0*(B*y + z)), 6.0*z + 6.0*A*x + 8.0*B*y + 8.0*z);
00232
00233 TESTER((A*x + B*y) + (A*y + B*x), (A + B)*x + (A+B)*y);
00234 TESTER((A*x + B*y) - (A*y + B*x), A*x - A*y + B*y - B*x);
00235
00236 TESTER((A*x + B*y) + 2.0*(A*y + B*x), A*(x + 2.0*y) + B*(2.0*x + y));
00237 TESTER((A*x + B*y) - 2.0*(A*y + B*x), A*(x - 2.0*y) + B*(y - 2.0*x));
00238
00239
00240 return pass;
00241 }
00242
00243
00244 template <class Scalar>
00245 inline bool LinearCombinationTester<Scalar>
00246 ::selfModifyingOpTests() const
00247 {
00248 bool pass = true;
00249
00250 RandomSparseMatrixBuilder<double> ABuilder(nLocalRows_, nLocalRows_,
00251 onProcDensity_, offProcDensity_,
00252 vecType_);
00253 RandomSparseMatrixBuilder<double> BBuilder(nLocalRows_, nLocalRows_,
00254 onProcDensity_, offProcDensity_,
00255 vecType_);
00256 RandomSparseMatrixBuilder<double> CBuilder(nLocalRows_, nLocalRows_,
00257 onProcDensity_, offProcDensity_,
00258 vecType_);
00259
00260
00261 LinearOperator<double> A = ABuilder.getOp();
00262 LinearOperator<double> B = BBuilder.getOp();
00263 LinearOperator<double> C = CBuilder.getOp();
00264
00265 VectorSpace<Scalar> vs = A.domain();
00266 A = identityOperator(vs);
00267 B = identityOperator(vs);
00268 C = identityOperator(vs);
00269
00270 Vector<Scalar> x = A.domain().createMember();
00271 Vector<Scalar> y = A.domain().createMember();
00272 Vector<Scalar> z = A.domain().createMember();
00273
00274 this->randomizeVec(x);
00275 this->randomizeVec(y);
00276 this->randomizeVec(z);
00277
00278 Vector<Scalar> a = x.copy();
00279 Vector<Scalar> b = y.copy();
00280 Vector<Scalar> c = z.copy();
00281
00282
00283 Out::os() << "starting linear combination tests" << std::endl;
00284
00285 x = 2.0*A*x;
00286 Scalar err = (x - 2.0*A*a).norm2();
00287 if (!this->checkTest(spec_, err, "x=2.0*A*x")) pass = false;
00288
00289 a = x.copy();
00290 x = x + y;
00291 err = (x - (a + y)).norm2();
00292 if (!this->checkTest(spec_, err, "x=x+y")) pass = false;
00293
00294 a = x.copy();
00295 x = y + x;
00296 err = (x - (y + a)).norm2();
00297 if (!this->checkTest(spec_, err, "x=y+x")) pass = false;
00298
00299 a = x.copy();
00300 x = A*x + B*y;
00301 err = (x - (A*a + B*y)).norm2();
00302 if (!this->checkTest(spec_, err, "x=A*x+B*y")) pass = false;
00303
00304 a = x.copy();
00305 x = B*y + A*x;
00306 err = (x - (B*y + A*a)).norm2();
00307 if (!this->checkTest(spec_, err, "x=B*y+A*x")) pass = false;
00308
00309 a = x.copy();
00310 x = A*x + (B*y + C*x);
00311 err = (x - (A*a + (B*y + C*a))).norm2();
00312 if (!this->checkTest(spec_, err, "x=A*x + (B*y + C*x)")) pass = false;
00313
00314 a = x.copy();
00315 x = (A*x + B*y) + C*x;
00316 err = (x - ((A*a + B*y) + C*a)).norm2();
00317 if (!this->checkTest(spec_, err, "x=(A*x + B*y) + C*x")) pass = false;
00318
00319
00320 Vector<Scalar> u;
00321 u = 2.0*A*B*x;
00322 err = (u - 2.0*A*B*x).norm2();
00323 if (!this->checkTest(spec_, err, "(empty) u=2*A*B*x")) pass = false;
00324
00325 u = 2.0*A*B*x;
00326 err = (u - 2.0*A*B*x).norm2();
00327 if (!this->checkTest(spec_, err, "(non-empty) u=2*A*B*x")) pass = false;
00328
00329
00330 Vector<Scalar> v;
00331 v = 2.0*x + 3.0*y;
00332 err = (v - (2.0*x + 3.0*y)).norm2();
00333 if (!this->checkTest(spec_, err, "(empty) v=2*x + 3*y")) pass = false;
00334
00335 v = 2.0*x + 3.0*y;
00336 err = (v - (2.0*x + 3.0*y)).norm2();
00337 if (!this->checkTest(spec_, err, "(non-empty) v=2*x + 3*y")) pass = false;
00338
00339
00340 Vector<Scalar> w;
00341 w = 2.0*x + 3.0*y + 5.0*z;
00342 err = (w - (2.0*x + 3.0*y + 5.0*z)).norm2();
00343 if (!this->checkTest(spec_, err, "(empty) w=2*x + 3*y + 5*z")) pass = false;
00344
00345 w = 2.0*x + 3.0*y + 5.0*z;
00346 err = (w - (2.0*x + 3.0*y + 5.0*z)).norm2();
00347 if (!this->checkTest(spec_, err, "(non-empty) w=2*x + 3*y + 5*z")) pass = false;
00348
00349
00350 Vector<Scalar> w2;
00351 w2 = 2.0*x + 3.0*y + 5.0*z + 7.0*u;
00352 err = (w2 - (2.0*x + 3.0*y + 5.0*z + 7.0*u)).norm2();
00353 if (!this->checkTest(spec_, err,
00354 "(empty) w2=2*x + 3*y + 5*z + 7*u")) pass = false;
00355
00356 w2 = 2.0*x + 3.0*y + 5.0*z + 7.0*u;
00357 err = (w2 - (2.0*x + 3.0*y + 5.0*z + 7.0*u)).norm2();
00358 if (!this->checkTest(spec_, err,
00359 "(non-empty) w2=2*x + 3*y + 5*z + 7*u")) pass = false;
00360
00361
00362 x = 2.0*x + 3.0*y + 5.0*z;
00363 err = (w - x).norm2();
00364 if (!this->checkTest(spec_, err, "x=2*x + 3*y + 5*z")) pass = false;
00365
00366
00367 return pass;
00368 }
00369
00370
00371 }
00372 #endif