Thyra Package Browser (Single Doxygen Collection) Version of the Day
EpetraExtAddTransformer_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 //
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 
00045 #include "Thyra_EpetraExtAddTransformer.hpp"
00046 #include "Thyra_DefaultScaledAdjointLinearOp.hpp"
00047 #include "Thyra_DefaultMultipliedLinearOp.hpp"
00048 #include "Thyra_DefaultAddedLinearOp.hpp"
00049 #include "Thyra_DefaultDiagonalLinearOp.hpp"
00050 #include "Thyra_VectorStdOps.hpp"
00051 #include "Thyra_TestingTools.hpp"
00052 #include "Thyra_LinearOpTester.hpp"
00053 #include "Thyra_EpetraThyraWrappers.hpp"
00054 #include "Thyra_EpetraLinearOp.hpp"
00055 #include "Thyra_get_Epetra_Operator.hpp"
00056 #include "Thyra_DefaultIdentityLinearOp.hpp"
00057 #include "EpetraExt_readEpetraLinearSystem.h"
00058 #include "Teuchos_GlobalMPISession.hpp"
00059 #include "Teuchos_VerboseObject.hpp"
00060 #include "Teuchos_XMLParameterListHelpers.hpp"
00061 #include "Teuchos_CommandLineProcessor.hpp"
00062 #include "Teuchos_StandardCatchMacros.hpp"
00063 
00064 #ifdef HAVE_MPI
00065 #  include "Epetra_MpiComm.h"
00066 #else
00067 #  include "Epetra_SerialComm.h"
00068 #endif
00069 
00070 #include "Teuchos_UnitTestHarness.hpp"
00071 
00072 
00073 namespace {
00074 
00075 
00076 using Teuchos::null;
00077 using Teuchos::RCP;
00078 using Thyra::EpetraExtAddTransformer;
00079 using Thyra::epetraExtAddTransformer;
00080 using Thyra::VectorBase;
00081 using Thyra::LinearOpBase;
00082 using Thyra::createMember;
00083 using Thyra::LinearOpTester;
00084 using Thyra::adjoint;
00085 using Thyra::multiply;
00086 using Thyra::diagonal;
00087 
00088 
00089 std::string matrixFile = "";
00090 std::string matrixFile2 = "";
00091 
00092 
00093 TEUCHOS_STATIC_SETUP()
00094 {
00095   Teuchos::UnitTestRepository::getCLP().setOption(
00096     "matrix-file", &matrixFile,
00097     "Defines the Epetra_CrsMatrix to read in."  );
00098   Teuchos::UnitTestRepository::getCLP().setOption(
00099     "matrix-file-2", &matrixFile2,
00100     "Defines the Epetra_CrsMatrix to read in."  );
00101 }
00102 
00103 const Teuchos::RCP<const Thyra::LinearOpBase<double> >
00104 buildAddOperator(int scenario,const Teuchos::RCP<const Thyra::LinearOpBase<double> > & A,
00105                                const Teuchos::RCP<const Thyra::LinearOpBase<double> > & B)
00106 {
00107    // build operators for the various addition/adjoint scenarios
00108    RCP<const Thyra::LinearOpBase<double> > M;
00109    
00110    switch(scenario) {
00111    case 0:
00112       M = Thyra::add(A,B,"A+B");
00113       break;
00114    case 1:
00115       M = Thyra::add(A,B,"A+adj(B)");
00116       break;
00117    case 2:
00118       M = Thyra::add(A,B,"adj(A)+B");
00119       break;
00120    case 3:
00121       M = Thyra::add(A,B,"adb(A)+adb(B)");
00122       break;
00123    default:
00124       TEUCHOS_ASSERT(false);
00125       break;
00126    }
00127 
00128    return M;
00129 }
00130 
00131 TEUCHOS_UNIT_TEST( EpetraExtAddTransformer, diag_mat_Add )
00132 {
00133   
00134   //
00135   // A) Read in problem matrices
00136   //
00137   
00138   out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\' ...\n";
00139     
00140 #ifdef HAVE_MPI
00141   Epetra_MpiComm comm(MPI_COMM_WORLD);
00142 #else
00143   Epetra_SerialComm comm;
00144 #endif
00145   RCP<Epetra_CrsMatrix> crsMat;
00146   EpetraExt::readEpetraLinearSystem( matrixFile, comm, &crsMat, NULL, NULL, NULL );
00147   
00148   //
00149   // B) Create the Thyra wrapped version
00150   //
00151   double scaleMat=3.7;
00152   //double scaleDiag=-2.9;
00153 
00154  
00155   const RCP<const Thyra::LinearOpBase<double> > A = Thyra::scale<double>(scaleMat,Thyra::epetraLinearOp(crsMat));
00156   RCP<VectorBase<double> > d = createMember(A->domain(), "d");
00157   V_S( d.ptr(), 3.0 ); // ToDo: Set ton != 1.0 and generalize
00158   const RCP<const Thyra::LinearOpBase<double> > B = diagonal(d);
00159 
00160   out << "\nA = " << *A;
00161   out << "\nB = " << *B;
00162 
00163   for(int scenario=0;scenario<4;scenario++) {
00164      //
00165      // C) Create implicit A+B operator
00166      //
00167    
00168      const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,A,B);
00169    
00170      //
00171      // D) Do the transformation
00172      //
00173    
00174      const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
00175    
00176      TEST_ASSERT(ApB_transformer != null);
00177    
00178      const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
00179      ApB_transformer->transform( *M, M_explicit.ptr() );
00180    
00181      out << "\nM_explicit = " << *M_explicit;
00182    
00183      //
00184      // E) Check the explicit operator
00185      //
00186    
00187      LinearOpTester<double> M_explicit_tester;
00188      M_explicit_tester.show_all_tests(true);;
00189    
00190      const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inOutArg(out) );
00191      if (!result) success = false;
00192   }
00193 
00194   for(int scenario=0;scenario<4;scenario++) {
00195      //
00196      // C) Create implicit A+B operator
00197      //
00198    
00199      const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,B,A);
00200    
00201      //
00202      // D) Do the transformation
00203      //
00204    
00205      const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
00206    
00207      TEST_ASSERT(ApB_transformer != null);
00208    
00209      const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
00210      ApB_transformer->transform( *M, M_explicit.ptr() );
00211    
00212      out << "\nM_explicit = " << *M_explicit;
00213    
00214      //
00215      // E) Check the explicit operator
00216      //
00217    
00218      LinearOpTester<double> M_explicit_tester;
00219      M_explicit_tester.show_all_tests(true);;
00220    
00221      const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inOutArg(out) );
00222      if (!result) success = false;
00223   }
00224 }
00225 
00226 TEUCHOS_UNIT_TEST( EpetraExtAddTransformer, id_mat_Add )
00227 {
00228   
00229   //
00230   // A) Read in problem matrices
00231   //
00232   
00233   out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\' ...\n";
00234     
00235 #ifdef HAVE_MPI
00236   Epetra_MpiComm comm(MPI_COMM_WORLD);
00237 #else
00238   Epetra_SerialComm comm;
00239 #endif
00240   RCP<Epetra_CrsMatrix> crsMat;
00241   EpetraExt::readEpetraLinearSystem( matrixFile, comm, &crsMat, NULL, NULL, NULL );
00242   
00243   //
00244   // B) Create the Thyra wrapped version
00245   //
00246   double scaleMat=3.7;
00247   //double scaleDiag=-2.9;
00248 
00249  
00250   const RCP<const Thyra::LinearOpBase<double> > A = Thyra::scale<double>(scaleMat,Thyra::epetraLinearOp(crsMat));
00251   const RCP<const Thyra::LinearOpBase<double> > B = identity(A->domain(),"id");
00252 
00253   out << "\nA = " << *A;
00254   out << "\nB = " << *B;
00255 
00256   for(int scenario=0;scenario<4;scenario++) {
00257      //
00258      // C) Create implicit A+B operator
00259      //
00260    
00261      const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,A,B);
00262    
00263      //
00264      // D) Do the transformation
00265      //
00266    
00267      const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
00268    
00269      TEST_ASSERT(ApB_transformer != null);
00270    
00271      const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
00272      ApB_transformer->transform( *M, M_explicit.ptr() );
00273    
00274      out << "\nM_explicit = " << *M_explicit;
00275    
00276      //
00277      // E) Check the explicit operator
00278      //
00279    
00280      LinearOpTester<double> M_explicit_tester;
00281      M_explicit_tester.show_all_tests(true);;
00282    
00283      const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inOutArg(out) );
00284      if (!result) success = false;
00285   }
00286 
00287   for(int scenario=0;scenario<4;scenario++) {
00288      //
00289      // C) Create implicit A+B operator
00290      //
00291    
00292      const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,B,A);
00293    
00294      //
00295      // D) Do the transformation
00296      //
00297    
00298      const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
00299    
00300      TEST_ASSERT(ApB_transformer != null);
00301    
00302      const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
00303      ApB_transformer->transform( *M, M_explicit.ptr() );
00304    
00305      out << "\nM_explicit = " << *M_explicit;
00306    
00307      //
00308      // E) Check the explicit operator
00309      //
00310    
00311      LinearOpTester<double> M_explicit_tester;
00312      M_explicit_tester.show_all_tests(true);;
00313    
00314      const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inOutArg(out) );
00315      if (!result) success = false;
00316   }
00317 }
00318 
00319 TEUCHOS_UNIT_TEST( EpetraExtAddTransformer, basic_Add )
00320 {
00321   
00322   //
00323   // A) Read in problem matrices
00324   //
00325   
00326   out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\' ...\n";
00327   out << "\nReading linear system in Epetra format from the file \'"<<matrixFile2<<"\' ...\n";
00328     
00329 #ifdef HAVE_MPI
00330   Epetra_MpiComm comm(MPI_COMM_WORLD);
00331 #else
00332   Epetra_SerialComm comm;
00333 #endif
00334   RCP<Epetra_CrsMatrix> epetra_A;
00335   RCP<Epetra_CrsMatrix> epetra_B;
00336   EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, NULL, NULL );
00337   EpetraExt::readEpetraLinearSystem( matrixFile2, comm, &epetra_B, NULL, NULL, NULL );
00338   
00339   //
00340   // B) Create the Thyra wrapped version
00341   //
00342   double scaleA=3.7;
00343   double scaleB=-2.9;
00344  
00345   const RCP<const Thyra::LinearOpBase<double> > A = Thyra::scale<double>(scaleA,Thyra::epetraLinearOp(epetra_B));
00346   const RCP<const Thyra::LinearOpBase<double> > B = Thyra::scale<double>(scaleB,Thyra::epetraLinearOp(epetra_B));
00347 
00348   out << "\nA = " << *A;
00349   out << "\nB = " << *B;
00350 
00351   for(int scenario=0;scenario<4;scenario++) {
00352      //
00353      // C) Create implicit A+B operator
00354      //
00355    
00356      const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,A,B);
00357    
00358      //
00359      // D) Do the transformation
00360      //
00361    
00362      const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
00363    
00364      TEST_ASSERT(ApB_transformer != null);
00365    
00366      const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
00367      ApB_transformer->transform( *M, M_explicit.ptr() );
00368    
00369      out << "\nM_explicit = " << *M_explicit;
00370    
00371      //
00372      // E) Check the explicit operator
00373      //
00374    
00375      LinearOpTester<double> M_explicit_tester;
00376      M_explicit_tester.show_all_tests(true);;
00377    
00378      const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inOutArg(out) );
00379      if (!result) success = false;
00380   }
00381 }
00382 
00383 TEUCHOS_UNIT_TEST( EpetraExtAddTransformer, mod_Add )
00384 {
00385   
00386   //
00387   // A) Read in problem matrices
00388   //
00389   
00390   out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\' ...\n";
00391   out << "\nReading linear system in Epetra format from the file \'"<<matrixFile2<<"\' ...\n";
00392     
00393 #ifdef HAVE_MPI
00394   Epetra_MpiComm comm(MPI_COMM_WORLD);
00395 #else
00396   Epetra_SerialComm comm;
00397 #endif
00398   RCP<Epetra_CrsMatrix> epetra_A;
00399   RCP<Epetra_CrsMatrix> epetra_B;
00400   EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, NULL, NULL );
00401   EpetraExt::readEpetraLinearSystem( matrixFile2, comm, &epetra_B, NULL, NULL, NULL );
00402   
00403   //
00404   // B) Create the Thyra wrapped version
00405   //
00406   double scaleA=3.7;
00407   double scaleB=-2.9;
00408  
00409   const RCP<const Thyra::LinearOpBase<double> > A = Thyra::scale<double>(scaleA,Thyra::epetraLinearOp(epetra_B));
00410   const RCP<const Thyra::LinearOpBase<double> > B = Thyra::scale<double>(scaleB,Thyra::epetraLinearOp(epetra_B));
00411 
00412   out << "\nA = " << *A;
00413   out << "\nB = " << *B;
00414 
00415   for(int scenario=0;scenario<4;scenario++) {
00416      //
00417      // C) Create implicit A+B operator
00418      //
00419    
00420      const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,A,B);
00421    
00422      //
00423      // D) Do the transformation
00424      //
00425    
00426      const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
00427    
00428      TEST_ASSERT(ApB_transformer != null);
00429    
00430      const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
00431      const RCP<LinearOpBase<double> > M_explicit_orig = M_explicit;
00432      ApB_transformer->transform( *M, M_explicit.ptr() );
00433 
00434      // do some violence to one of the operators
00435      double * view; int numEntries;
00436      epetra_B->Scale(3.2);
00437      TEUCHOS_ASSERT(epetra_B->ExtractMyRowView(3,numEntries,view)==0);
00438      for(int i=0;i<numEntries;i++) view[i] += view[i]*double(i+1);
00439 
00440      // compute multiply again
00441      ApB_transformer->transform( *M, M_explicit.ptr() );
00442 
00443      out << "\nM_explicit = " << *M_explicit;
00444 
00445      TEUCHOS_ASSERT(Thyra::get_Epetra_Operator(*M_explicit)
00446                   ==Thyra::get_Epetra_Operator(*M_explicit_orig));
00447    
00448      //
00449      // E) Check the explicit operator
00450      //
00451    
00452      LinearOpTester<double> M_explicit_tester;
00453      M_explicit_tester.show_all_tests(true);;
00454    
00455      const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inoutArg(out) );
00456      if (!result) success = false;
00457   }
00458 }
00459 
00460 } // end namespace
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines