Thyra Package Browser (Single Doxygen Collection) Version of the Day
EpetraLinearOp_UnitTests.cpp
Go to the documentation of this file.
00001 /*
00002 // @HEADER
00003 // ***********************************************************************
00004 // 
00005 //    Thyra: Interfaces and Support for Abstract Numerical Algorithms
00006 //                 Copyright (2004) Sandia Corporation
00007 // 
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 // 
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 1//
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov) 
00039 // 
00040 // ***********************************************************************
00041 // @HEADER
00042 */
00043 
00044 #include "Thyra_EpetraLinearOp.hpp"
00045 #include "Thyra_ScaledLinearOpBase.hpp"
00046 #include "Thyra_RowStatLinearOpBase.hpp"
00047 #include "Thyra_LinearOpTester.hpp"
00048 #include "Thyra_DefaultBlockedLinearOp.hpp"
00049 #include "Thyra_MultiVectorBase.hpp"
00050 #include "Thyra_MultiVectorStdOps.hpp"
00051 #include "EpetraThyraAdaptersTestHelpers.hpp"
00052 
00053 #include "Teuchos_UnitTestHarness.hpp"
00054 
00055 
00056 namespace {
00057 
00058 
00059 bool g_dumpAll = false;
00060 const int g_dim = 4;
00061 
00062 
00063 TEUCHOS_STATIC_SETUP()
00064 {
00065   Teuchos::UnitTestRepository::getCLP().setOption(
00066     "dump-all", "no-dump-all", &g_dumpAll,
00067     "Dump lots of data" );
00068 }
00069 
00070 
00071 } // namespace
00072 
00073 
00074 namespace Thyra {
00075 
00076 
00077 //
00078 // Unit Tests
00079 //
00080 
00081 
00082 TEUCHOS_UNIT_TEST( EpetraLinearOp, ScaledLinearOpBase )
00083 {
00084   using Teuchos::null;
00085   using Teuchos::inOutArg;
00086   using Teuchos::updateSuccess;
00087   using Teuchos::rcp_dynamic_cast;
00088   typedef ScalarTraits<double> ST;
00089 
00090   // Set up an EpetraLinearOp
00091 
00092   const RCP<const Epetra_Comm> comm = getEpetraComm();
00093   const int numLocalRows = g_localDim;
00094   const int numRows = numLocalRows * comm->NumProc();
00095   const int numCols = numLocalRows / 2;
00096   const RCP<Epetra_CrsMatrix> epetraCrsM = getEpetraMatrix(numRows, numCols);
00097   const RCP<LinearOpBase<double> > epetraOp = nonconstEpetraLinearOp(epetraCrsM);
00098 
00099   const RCP<ScaledLinearOpBase<double> > scaledOp =
00100     rcp_dynamic_cast<ScaledLinearOpBase<double> >(epetraOp, true);
00101 
00102   // Get the original mat-vec
00103 
00104   const double two = 2.0;
00105 
00106   const RCP<VectorBase<double> > rhs_vec =
00107     createMember<double>(epetraOp->domain());
00108   assign<double>(rhs_vec.ptr(), two);
00109 
00110   const RCP<VectorBase<double> > lhs_orig_vec =
00111     createMember<double>(epetraOp->range());
00112 
00113   apply<double>(*epetraOp, NOTRANS, *rhs_vec, lhs_orig_vec.ptr());
00114 
00115   if (g_dumpAll) {
00116     out << "epetraOp = " << *epetraOp;
00117     out << "rhs_vec = " << *rhs_vec;
00118     out << "lhs_orig_vec = " << *lhs_orig_vec;
00119   }
00120 
00121   // Scale the op from the left (row scaling)
00122 
00123   const double three = 3.0;
00124 
00125   const RCP<VectorBase<double> > row_scaling =
00126     createMember<double>(epetraOp->range());
00127   assign<double>(row_scaling.ptr(), three);
00128 
00129   scaledOp->scaleLeft(*row_scaling);
00130 
00131   if (g_dumpAll) {
00132     out << "row_scaling = " << *row_scaling;
00133     out << "epetraOp left scaled = " << *epetraOp;
00134   }
00135 
00136   // Test that resulting left scaling
00137 
00138   const RCP<VectorBase<double> > lhs_left_scaled_vec =
00139     createMember<double>(epetraOp->range());
00140 
00141   apply<double>(*epetraOp, NOTRANS, *rhs_vec, lhs_left_scaled_vec.ptr());
00142 
00143   if (g_dumpAll) {
00144     out << "lhs_left_scaled_vec = " << *lhs_left_scaled_vec;
00145   }
00146 
00147   TEST_FLOATING_EQUALITY(
00148     three * sum<double>(*lhs_orig_vec),
00149     sum<double>(*lhs_left_scaled_vec),
00150     as<double>(10.0 * ST::eps())
00151     );
00152 
00153   // Left scale the matrix back
00154 
00155   const RCP<VectorBase<double> > inv_row_scaling =
00156     createMember<double>(epetraOp->range());
00157   reciprocal<double>(*row_scaling, inv_row_scaling.ptr());
00158 
00159   scaledOp->scaleLeft(*inv_row_scaling);
00160 
00161   if (g_dumpAll) {
00162     out << "inv_row_scaling = " << *row_scaling;
00163     out << "epetraOp left scaled back to orig = " << *epetraOp;
00164   }
00165 
00166   const RCP<VectorBase<double> > lhs_orig2_vec =
00167     createMember<double>(epetraOp->range());
00168 
00169   apply<double>(*epetraOp, NOTRANS, *rhs_vec, lhs_orig2_vec.ptr());
00170 
00171   if (g_dumpAll) {
00172     out << "lhs_orig2_vec = " << *lhs_orig2_vec;
00173   }
00174 
00175   TEST_FLOATING_EQUALITY(
00176     sum<double>(*lhs_orig_vec),
00177     sum<double>(*lhs_orig2_vec),
00178     as<double>(10.0 * ST::eps())
00179     );
00180   // NOTE: Above, it would ask for exact binary match except if one uses
00181   // threading it will not match exactly!
00182 
00183   // Scale the op from the right (col scaling)
00184 
00185   const double four = 4.0;
00186 
00187   const RCP<VectorBase<double> > col_scaling =
00188     createMember<double>(epetraOp->domain());
00189   assign<double>(col_scaling.ptr(), four);
00190 
00191   scaledOp->scaleRight(*col_scaling);
00192 
00193   if (g_dumpAll) {
00194     out << "col_scaling = " << *col_scaling;
00195     out << "epetraOp right scaled = " << *epetraOp;
00196   }
00197 
00198   // Test that resulting right scaling
00199 
00200   const RCP<VectorBase<double> > lhs_right_scaled_vec =
00201     createMember<double>(epetraOp->range());
00202 
00203   apply<double>(*epetraOp, NOTRANS, *rhs_vec, lhs_right_scaled_vec.ptr());
00204 
00205   if (g_dumpAll) {
00206     out << "lhs_right_scaled_vec = " << *lhs_right_scaled_vec;
00207   }
00208 
00209   TEST_FLOATING_EQUALITY(
00210     four * sum<double>(*lhs_orig_vec),
00211     sum<double>(*lhs_right_scaled_vec),
00212     as<double>(10.0 * ST::eps())
00213     );
00214 
00215   // Right scale the matrix back
00216 
00217   const RCP<VectorBase<double> > inv_col_scaling =
00218     createMember<double>(epetraOp->domain());
00219   reciprocal<double>(*col_scaling, inv_col_scaling.ptr());
00220 
00221   scaledOp->scaleRight(*inv_col_scaling);
00222 
00223   if (g_dumpAll) {
00224     out << "inv_col_scaling = " << *col_scaling;
00225     out << "epetraOp right scaled back to orig = " << *epetraOp;
00226   }
00227 
00228   const RCP<VectorBase<double> > lhs_orig3_vec =
00229     createMember<double>(epetraOp->range());
00230 
00231   apply<double>(*epetraOp, NOTRANS, *rhs_vec, lhs_orig3_vec.ptr());
00232 
00233   if (g_dumpAll) {
00234     out << "lhs_orig3_vec = " << *lhs_orig3_vec;
00235   }
00236 
00237   TEST_FLOATING_EQUALITY(
00238     sum<double>(*lhs_orig_vec),
00239     sum<double>(*lhs_orig3_vec),
00240     as<double>(10.0 * ST::eps())
00241     );
00242   // NOTE: Above, it would ask for exact binary match except if one uses
00243   // threading it will not match exactly!
00244   
00245 }
00246 
00247 
00248 TEUCHOS_UNIT_TEST( EpetraLinearOp, RowStatLinearOpBase )
00249 {
00250   using Teuchos::null;
00251   using Teuchos::inOutArg;
00252   using Teuchos::updateSuccess;
00253   using Teuchos::rcp_dynamic_cast;
00254   typedef ScalarTraits<double> ST;
00255 
00256   // Set up the EpetraLinearOp
00257 
00258   const RCP<const Epetra_Comm> comm = getEpetraComm();
00259   const int numLocalRows = g_localDim;
00260   const int numRows = numLocalRows * comm->NumProc();
00261   const int numCols = numLocalRows / 2;
00262   const RCP<Epetra_CrsMatrix> epetraCrsM = getEpetraMatrix(numRows, numCols);
00263   const double two = 2.0;
00264   epetraCrsM->PutScalar(two);
00265   const RCP<LinearOpBase<double> > epetraOp = nonconstEpetraLinearOp(epetraCrsM);
00266 
00267   const RCP<RowStatLinearOpBase<double> > rowStatOp =
00268     rcp_dynamic_cast<RowStatLinearOpBase<double> >(epetraOp, true);
00269 
00270   if (g_dumpAll) {
00271     out << "epetraOp = " << *epetraOp;
00272   }
00273 
00274   // Get the inverse row sums
00275 
00276   const RCP<VectorBase<double> > inv_row_sums =
00277     createMember<double>(epetraOp->range());
00278 
00279   rowStatOp->getRowStat(RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM,
00280     inv_row_sums.ptr());
00281 
00282   if (g_dumpAll) {
00283     out << "inv_row_sums = " << *inv_row_sums;
00284   }
00285 
00286   TEST_FLOATING_EQUALITY(
00287     sum<double>(*inv_row_sums),
00288     as<double>((1.0 / (two * numCols)) * numRows),
00289     as<double>(10.0 * ST::eps())
00290     );
00291 
00292 }
00293 
00294 
00295 TEUCHOS_UNIT_TEST( EpetraLinearOp, rectangular )
00296 {
00297   using Teuchos::null;
00298   using Teuchos::inOutArg;
00299   using Teuchos::updateSuccess;
00300 
00301   const RCP<const Epetra_Comm> comm = getEpetraComm();
00302 
00303   const int numLocalRows = g_localDim;
00304   const int numRows = numLocalRows * comm->NumProc();
00305   const int numCols = numLocalRows / 2;
00306 
00307   const RCP<Epetra_CrsMatrix> epetraCrsM = getEpetraMatrix(numRows, numCols);
00308 
00309   const RCP<const LinearOpBase<double> > epetraOp = epetraLinearOp(epetraCrsM);
00310 
00311   LinearOpTester<double> linearOpTester;
00312   updateSuccess(linearOpTester.check(*epetraOp, inOutArg(out)), success);
00313    
00314 }
00315 
00316 
00317 TEUCHOS_UNIT_TEST( EpetraLinearOp, blocked_op )
00318 {
00319 
00320   using Teuchos::describe;
00321   
00322   // build sub operators
00323   RCP<const LinearOpBase<double> > A00 =
00324     epetraLinearOp(getEpetraMatrix(4,4,0));
00325   RCP<const LinearOpBase<double> > A01 =
00326     epetraLinearOp(getEpetraMatrix(4,3,1));
00327   RCP<const LinearOpBase<double> > A02 =
00328     epetraLinearOp(getEpetraMatrix(4,2,2));
00329   RCP<const LinearOpBase<double> > A10 =
00330     epetraLinearOp(getEpetraMatrix(3,4,3));
00331   RCP<const LinearOpBase<double> > A11 =
00332     epetraLinearOp(getEpetraMatrix(3,3,4));
00333   RCP<const LinearOpBase<double> > A12 =
00334     epetraLinearOp(getEpetraMatrix(3,2,5));
00335   RCP<const LinearOpBase<double> > A20 =
00336     epetraLinearOp(getEpetraMatrix(2,4,6));
00337   RCP<const LinearOpBase<double> > A21 =
00338     epetraLinearOp(getEpetraMatrix(2,3,8));
00339   RCP<const LinearOpBase<double> > A22 =
00340     epetraLinearOp(getEpetraMatrix(2,2,8));
00341   
00342   out << "Sub operators built" << std::endl;
00343 
00344   {
00345      // build composite operator
00346      RCP<const LinearOpBase<double> > A = block2x2<double>(
00347        block2x2<double>(A00, A01, A10, A11),   block2x1<double>(A02,A12),
00348        block1x2<double>(A20, A21),             A22
00349        );
00350    
00351      out << "First composite operator built" << std::endl;
00352      
00353      // build vectors for use in apply
00354      RCP<MultiVectorBase<double> > x = createMembers<double>(A->domain(), 3);
00355      RCP<MultiVectorBase<double> > y = createMembers<double>(A->range(), 3);
00356      
00357      randomize(-1.0, 1.0, x.ptr());
00358    
00359      out << "A = \n" << describe(*A, Teuchos::VERB_HIGH) << std::endl;
00360      out << "x = \n" << describe(*x, Teuchos::VERB_HIGH) << std::endl;
00361      out << "y = \n" << describe(*y, Teuchos::VERB_HIGH) << std::endl;
00362      
00363      // perform a matrix vector multiply
00364      apply(*A, NOTRANS, *x, y.ptr());
00365 
00366      out << "First composite operator completed" << std::endl;
00367   }
00368 
00369   {
00370      RCP<const LinearOpBase<double> > A = block2x2<double>(
00371        A11,                          block1x2<double>(A10, A12),
00372        block2x1<double>(A01, A21),   block2x2<double>(A00, A02, A20, A22)
00373        );
00374      
00375      out << "Second composite operator built" << std::endl;
00376      
00377      // build vectors for use in apply
00378      RCP<MultiVectorBase<double> > x = createMembers<double>(A->domain(), 3);
00379      RCP<MultiVectorBase<double> > y = createMembers<double>(A->range(), 3);
00380      
00381      randomize(-1.0, 1.0, x.ptr());
00382    
00383      out << "A = \n" << describe(*A, Teuchos::VERB_HIGH) << std::endl;
00384      out << "x = \n" << describe(*x, Teuchos::VERB_HIGH) << std::endl;
00385      out << "y = \n" << describe(*y, Teuchos::VERB_HIGH) << std::endl;
00386      
00387      // perform a matrix vector multiply
00388      apply(*A, NOTRANS, *x, y.ptr());
00389 
00390      out << "Second composite operator completed" << std::endl;
00391   }
00392 
00393   out << "Test complete" << std::endl;
00394 
00395 }
00396 
00397 
00398 } // namespace Thyra
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines