00001
00002
00003
00004
00005
00006 #ifndef PLAYA_MATRIXMATRIXTESTER_HPP
00007 #define PLAYA_MATRIXMATRIXTESTER_HPP
00008
00009 #include "PlayaLinearOperatorDecl.hpp"
00010 #include "PlayaEpetraMatrixMatrixProduct.hpp"
00011 #include "PlayaEpetraMatrixMatrixSum.hpp"
00012 #include "PlayaEpetraMatrixOps.hpp"
00013 #include "PlayaEpetraMatrix.hpp"
00014 #include "Teuchos_ScalarTraits.hpp"
00015 #include "PlayaLinearCombinationImpl.hpp"
00016
00017 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00018 #include "PlayaSimpleComposedOpImpl.hpp"
00019 #include "PlayaSimpleScaledOpImpl.hpp"
00020 #include "PlayaSimpleAddedOpImpl.hpp"
00021 #include "PlayaSimpleDiagonalOpImpl.hpp"
00022 #include "PlayaSimpleTransposedOpImpl.hpp"
00023 #include "PlayaRandomSparseMatrixBuilderImpl.hpp"
00024 #endif
00025
00026 using namespace Playa;
00027 using namespace Teuchos;
00028
00029
00030
00031 namespace Playa
00032 {
00033
00034
00035 template <class Scalar>
00036 class MatrixMatrixTester : public TesterBase<Scalar>
00037 {
00038 public:
00039
00040 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
00041
00042
00043 MatrixMatrixTester(const LinearOperator<Scalar>& A,
00044 const LinearOperator<Scalar>& B,
00045 const TestSpecifier<Scalar>& sumSpec,
00046 const TestSpecifier<Scalar>& prodSpec,
00047 const TestSpecifier<Scalar>& diagSpec,
00048 const TestSpecifier<Scalar>& diagLeftProdSpec,
00049 const TestSpecifier<Scalar>& diagRightProdSpec);
00050
00051
00052 bool runAllTests() const ;
00053
00054
00055 bool sumTest() const ;
00056
00057
00058 bool prodTest() const ;
00059
00060
00061 bool diagTest() const ;
00062
00063
00064 bool diagLeftProdTest() const ;
00065
00066
00067 bool diagRightProdTest() const ;
00068
00069
00070 private:
00071
00072 LinearOperator<Scalar> A_;
00073
00074 LinearOperator<Scalar> B_;
00075
00076 TestSpecifier<Scalar> sumSpec_;
00077
00078 TestSpecifier<Scalar> prodSpec_;
00079
00080 TestSpecifier<Scalar> diagSpec_;
00081
00082 TestSpecifier<Scalar> diagLeftProdSpec_;
00083
00084 TestSpecifier<Scalar> diagRightProdSpec_;
00085
00086 };
00087
00088 template <class Scalar>
00089 inline MatrixMatrixTester<Scalar>
00090 ::MatrixMatrixTester(const LinearOperator<Scalar>& A,
00091 const LinearOperator<Scalar>& B,
00092 const TestSpecifier<Scalar>& sumSpec,
00093 const TestSpecifier<Scalar>& prodSpec,
00094 const TestSpecifier<Scalar>& diagSpec,
00095 const TestSpecifier<Scalar>& diagRightProdSpec,
00096 const TestSpecifier<Scalar>& diagLeftProdSpec)
00097 : TesterBase<Scalar>(),
00098 A_(A),
00099 B_(B),
00100 sumSpec_(sumSpec),
00101 prodSpec_(prodSpec),
00102 diagSpec_(diagSpec),
00103 diagLeftProdSpec_(diagLeftProdSpec),
00104 diagRightProdSpec_(diagRightProdSpec)
00105 {;}
00106
00107 template <class Scalar>
00108 inline bool MatrixMatrixTester<Scalar>
00109 ::runAllTests() const
00110 {
00111 bool pass = true;
00112
00113 pass = this->sumTest() && pass;
00114 pass = this->prodTest() && pass;
00115 pass = this->diagTest() && pass;
00116 pass = this->diagLeftProdTest() && pass;
00117 pass = this->diagRightProdTest() && pass;
00118
00119 return pass;
00120 }
00121
00122 template <class Scalar>
00123 inline bool MatrixMatrixTester<Scalar>
00124 ::diagTest() const
00125 {
00126 if (diagSpec_.doTest())
00127 {
00128 Vector<Scalar> x = B_.domain().createMember();
00129 Vector<Scalar> d = B_.domain().createMember();
00130 this->randomizeVec(x);
00131 this->randomizeVec(d);
00132 LinearOperator<Scalar> D0 = diagonalOperator(d);
00133 LinearOperator<Scalar> D = makeEpetraDiagonalMatrix(d);
00134 Vector<Scalar> d1 = getEpetraDiagonal(D);
00135
00136 Out::root() << "computing implicit product y1 = D*x..." << std::endl;
00137 Vector<Scalar> y1 = D0*x;
00138 Out::root() << "computing explicit product y2 = D*x..." << std::endl;
00139 Vector<Scalar> y2 = D*x;
00140
00141 ScalarMag err = (y1 - y2).norm2();
00142
00143 Out::root() << "|y1-y2| = " << err << std::endl;
00144
00145 Out::root() << "comparing recovered and original diagonals" << std::endl;
00146 ScalarMag err2 = (d - d1).norm2();
00147 Out::root() << "|d1-d2| = " << err2 << std::endl;
00148
00149 return this->checkTest(prodSpec_, err+err2, "matrix-matrix multiply");
00150
00151 }
00152 Out::root() << "skipping matrix-matrix multiply test..." << std::endl;
00153 return true;
00154 }
00155
00156
00157
00158 template <class Scalar>
00159 inline bool MatrixMatrixTester<Scalar>
00160 ::sumTest() const
00161 {
00162 if (sumSpec_.doTest())
00163 {
00164
00165
00166 if (A_.range() != B_.range() || A_.domain() != B_.domain())
00167 {
00168 Out::root() << "skipping sum on incompatible matrices" << std::endl;
00169 return true;
00170 }
00171
00172 Out::root() << "running matrix-matrix multiply test..." << std::endl;
00173 LinearOperator<Scalar> implicitAdd = A_ + B_;
00174 LinearOperator<Scalar> explicitAdd = epetraMatrixMatrixSum(A_, B_);
00175
00176 Vector<Scalar> x = B_.domain().createMember();
00177 this->randomizeVec(x);
00178 Out::root() << "computing implicit sum y1 = (A+B)*x..." << std::endl;
00179 Vector<Scalar> y1 = implicitAdd*x;
00180 Out::root() << "computing explicit sum y2 = (A+B)*x..." << std::endl;
00181 Vector<Scalar> y2 = explicitAdd*x;
00182
00183 ScalarMag err = (y1 - y2).norm2();
00184
00185 Out::root() << "|y1-y2| = " << err << std::endl;
00186 return this->checkTest(prodSpec_, err, "matrix-matrix multiply");
00187
00188 }
00189 Out::root() << "skipping matrix-matrix multiply test..." << std::endl;
00190 return true;
00191 }
00192
00193
00194 template <class Scalar>
00195 inline bool MatrixMatrixTester<Scalar>
00196 ::prodTest() const
00197 {
00198 if (prodSpec_.doTest())
00199 {
00200 Out::root() << "running matrix-matrix multiply test..." << std::endl;
00201 LinearOperator<Scalar> composed = A_ * B_;
00202 LinearOperator<Scalar> multiplied = epetraMatrixMatrixProduct(A_, B_);
00203
00204 Vector<Scalar> x = B_.domain().createMember();
00205 this->randomizeVec(x);
00206 Out::root() << "computing implicit product y1 = (A*B)*x..." << std::endl;
00207 Vector<Scalar> y1 = composed*x;
00208 Out::root() << "computing explicit product y2 = (A*B)*x..." << std::endl;
00209 Vector<Scalar> y2 = multiplied*x;
00210
00211 ScalarMag err = (y1 - y2).norm2();
00212
00213 Out::root() << "|y1-y2| = " << err << std::endl;
00214 return this->checkTest(prodSpec_, err, "matrix-matrix multiply");
00215 }
00216 Out::root() << "skipping matrix-matrix multiply test..." << std::endl;
00217 return true;
00218 }
00219
00220
00221 template <class Scalar>
00222 inline bool MatrixMatrixTester<Scalar>
00223 ::diagLeftProdTest() const
00224 {
00225 if (diagLeftProdSpec_.doTest())
00226 {
00227 Out::root() << "running diagonal*matrix multiplication test..." << std::endl;
00228
00229 Vector<Scalar> x = A_.domain().createMember();
00230 this->randomizeVec(x);
00231
00232 Vector<Scalar> d = A_.range().createMember();
00233 this->randomizeVec(d);
00234
00235 LinearOperator<Scalar> D = diagonalOperator(d);
00236 LinearOperator<Scalar> DA = epetraLeftScale(d, A_);
00237
00238 Out::root() << "computing implicit y1 = D*A*x..." << std::endl;
00239 Vector<Scalar> y1 = D*A_*x;
00240 Out::root() << "computing explicit y2 = D*A*x..." << std::endl;
00241 Vector<Scalar> y2 = DA*x;
00242
00243 ScalarMag err = (y1 - y2).norm2();
00244
00245 Out::root() << "|y1-y2| = " << err << std::endl;
00246
00247 return this->checkTest(diagLeftProdSpec_, err, "diagonal*matrix multiplication");
00248 }
00249 Out::root() << "skipping diagonal matrix-matrix test..." << std::endl;
00250 return true;
00251 }
00252
00253
00254 template <class Scalar>
00255 inline bool MatrixMatrixTester<Scalar>
00256 ::diagRightProdTest() const
00257 {
00258 if (diagRightProdSpec_.doTest())
00259 {
00260 Out::root() << "running diagonal*matrix multiplication test..." << std::endl;
00261
00262 Vector<Scalar> x = A_.domain().createMember();
00263 this->randomizeVec(x);
00264
00265 Vector<Scalar> d = A_.domain().createMember();
00266 this->randomizeVec(d);
00267
00268 LinearOperator<Scalar> D = diagonalOperator(d);
00269 LinearOperator<Scalar> AD = epetraRightScale(A_, d);
00270
00271 Out::root() << "computing implicit y1 = A*D*x..." << std::endl;
00272 Vector<Scalar> y1 = A_*D*x;
00273 Out::root() << "computing explicit y2 = A*D*x..." << std::endl;
00274 Vector<Scalar> y2 = AD*x;
00275
00276 ScalarMag err = (y1 - y2).norm2();
00277
00278 Out::root() << "|y1-y2| = " << err << std::endl;
00279
00280 return this->checkTest(diagLeftProdSpec_, err, "matrix*diagonal multiplication");
00281 }
00282 Out::root() << "skipping diagonal matrix-matrix test..." << std::endl;
00283 return true;
00284 }
00285
00286
00287
00288 }
00289 #endif