Thyra Package Browser (Single Doxygen Collection) Version of the Day
EpetraExtAddTransformer_UnitTests.cpp
Go to the documentation of this file.
00001 
00002 #include "Thyra_EpetraExtAddTransformer.hpp"
00003 #include "Thyra_DefaultScaledAdjointLinearOp.hpp"
00004 #include "Thyra_DefaultMultipliedLinearOp.hpp"
00005 #include "Thyra_DefaultAddedLinearOp.hpp"
00006 #include "Thyra_DefaultDiagonalLinearOp.hpp"
00007 #include "Thyra_VectorStdOps.hpp"
00008 #include "Thyra_TestingTools.hpp"
00009 #include "Thyra_LinearOpTester.hpp"
00010 #include "Thyra_EpetraThyraWrappers.hpp"
00011 #include "Thyra_EpetraLinearOp.hpp"
00012 #include "Thyra_get_Epetra_Operator.hpp"
00013 #include "EpetraExt_readEpetraLinearSystem.h"
00014 #include "Teuchos_GlobalMPISession.hpp"
00015 #include "Teuchos_VerboseObject.hpp"
00016 #include "Teuchos_XMLParameterListHelpers.hpp"
00017 #include "Teuchos_CommandLineProcessor.hpp"
00018 #include "Teuchos_StandardCatchMacros.hpp"
00019 
00020 #ifdef HAVE_MPI
00021 #  include "Epetra_MpiComm.h"
00022 #else
00023 #  include "Epetra_SerialComm.h"
00024 #endif
00025 
00026 #include "Teuchos_UnitTestHarness.hpp"
00027 
00028 
00029 namespace {
00030 
00031 
00032 using Teuchos::null;
00033 using Teuchos::RCP;
00034 using Thyra::EpetraExtAddTransformer;
00035 using Thyra::epetraExtAddTransformer;
00036 using Thyra::VectorBase;
00037 using Thyra::LinearOpBase;
00038 using Thyra::createMember;
00039 using Thyra::LinearOpTester;
00040 using Thyra::adjoint;
00041 using Thyra::multiply;
00042 using Thyra::diagonal;
00043 
00044 
00045 std::string matrixFile = "";
00046 std::string matrixFile2 = "";
00047 
00048 
00049 TEUCHOS_STATIC_SETUP()
00050 {
00051   Teuchos::UnitTestRepository::getCLP().setOption(
00052     "matrix-file", &matrixFile,
00053     "Defines the Epetra_CrsMatrix to read in."  );
00054   Teuchos::UnitTestRepository::getCLP().setOption(
00055     "matrix-file-2", &matrixFile2,
00056     "Defines the Epetra_CrsMatrix to read in."  );
00057 }
00058 
00059 const Teuchos::RCP<const Thyra::LinearOpBase<double> >
00060 buildAddOperator(int scenario,const Teuchos::RCP<const Thyra::LinearOpBase<double> > & A,
00061                                const Teuchos::RCP<const Thyra::LinearOpBase<double> > & B)
00062 {
00063    // build operators for the various addition/adjoint scenarios
00064    RCP<const Thyra::LinearOpBase<double> > M;
00065    
00066    switch(scenario) {
00067    case 0:
00068       M = Thyra::add(A,B,"A+B");
00069       break;
00070    case 1:
00071       M = Thyra::add(A,B,"A+adj(B)");
00072       break;
00073    case 2:
00074       M = Thyra::add(A,B,"adj(A)+B");
00075       break;
00076    case 3:
00077       M = Thyra::add(A,B,"adb(A)+adb(B)");
00078       break;
00079    default:
00080       TEUCHOS_ASSERT(false);
00081       break;
00082    }
00083 
00084    return M;
00085 }
00086 
00087 TEUCHOS_UNIT_TEST( EpetraExtAddTransformer, basic_Add )
00088 {
00089   
00090   //
00091   // A) Read in problem matrices
00092   //
00093   
00094   out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\' ...\n";
00095   out << "\nReading linear system in Epetra format from the file \'"<<matrixFile2<<"\' ...\n";
00096     
00097 #ifdef HAVE_MPI
00098   Epetra_MpiComm comm(MPI_COMM_WORLD);
00099 #else
00100   Epetra_SerialComm comm;
00101 #endif
00102   RCP<Epetra_CrsMatrix> epetra_A;
00103   RCP<Epetra_CrsMatrix> epetra_B;
00104   EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, NULL, NULL );
00105   EpetraExt::readEpetraLinearSystem( matrixFile2, comm, &epetra_B, NULL, NULL, NULL );
00106   
00107   //
00108   // B) Create the Thyra wrapped version
00109   //
00110   double scaleA=3.7;
00111   double scaleB=-2.9;
00112  
00113   const RCP<const Thyra::LinearOpBase<double> > A = Thyra::scale<double>(scaleA,Thyra::epetraLinearOp(epetra_B));
00114   const RCP<const Thyra::LinearOpBase<double> > B = Thyra::scale<double>(scaleB,Thyra::epetraLinearOp(epetra_B));
00115 
00116   out << "\nA = " << *A;
00117   out << "\nB = " << *B;
00118 
00119   for(int scenario=0;scenario<4;scenario++) {
00120      //
00121      // C) Create implicit A+B operator
00122      //
00123    
00124      const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,A,B);
00125    
00126      //
00127      // D) Do the transformation
00128      //
00129    
00130      const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
00131    
00132      TEST_ASSERT(ApB_transformer != null);
00133    
00134      const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
00135      ApB_transformer->transform( *M, M_explicit.ptr() );
00136    
00137      out << "\nM_explicit = " << *M_explicit;
00138    
00139      //
00140      // E) Check the explicit operator
00141      //
00142    
00143      LinearOpTester<double> M_explicit_tester;
00144      M_explicit_tester.show_all_tests(true);;
00145    
00146      const bool result = M_explicit_tester.compare( *M, *M_explicit, &out );
00147      if (!result) success = false;
00148   }
00149 }
00150 
00151 TEUCHOS_UNIT_TEST( EpetraExtAddTransformer, mod_Add )
00152 {
00153   
00154   //
00155   // A) Read in problem matrices
00156   //
00157   
00158   out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\' ...\n";
00159   out << "\nReading linear system in Epetra format from the file \'"<<matrixFile2<<"\' ...\n";
00160     
00161 #ifdef HAVE_MPI
00162   Epetra_MpiComm comm(MPI_COMM_WORLD);
00163 #else
00164   Epetra_SerialComm comm;
00165 #endif
00166   RCP<Epetra_CrsMatrix> epetra_A;
00167   RCP<Epetra_CrsMatrix> epetra_B;
00168   EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, NULL, NULL );
00169   EpetraExt::readEpetraLinearSystem( matrixFile2, comm, &epetra_B, NULL, NULL, NULL );
00170   
00171   //
00172   // B) Create the Thyra wrapped version
00173   //
00174   double scaleA=3.7;
00175   double scaleB=-2.9;
00176  
00177   const RCP<const Thyra::LinearOpBase<double> > A = Thyra::scale<double>(scaleA,Thyra::epetraLinearOp(epetra_B));
00178   const RCP<const Thyra::LinearOpBase<double> > B = Thyra::scale<double>(scaleB,Thyra::epetraLinearOp(epetra_B));
00179 
00180   out << "\nA = " << *A;
00181   out << "\nB = " << *B;
00182 
00183   for(int scenario=0;scenario<4;scenario++) {
00184      //
00185      // C) Create implicit A+B operator
00186      //
00187    
00188      const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,A,B);
00189    
00190      //
00191      // D) Do the transformation
00192      //
00193    
00194      const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
00195    
00196      TEST_ASSERT(ApB_transformer != null);
00197    
00198      const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
00199      const RCP<LinearOpBase<double> > M_explicit_orig = M_explicit;
00200      ApB_transformer->transform( *M, M_explicit.ptr() );
00201 
00202      // do some violence to one of the operators
00203      double * view; int numEntries;
00204      epetra_B->Scale(3.2);
00205      TEUCHOS_ASSERT(epetra_B->ExtractMyRowView(3,numEntries,view)==0);
00206      for(int i=0;i<numEntries;i++) view[i] += view[i]*double(i+1);
00207 
00208      // compute multiply again
00209      ApB_transformer->transform( *M, M_explicit.ptr() );
00210 
00211      out << "\nM_explicit = " << *M_explicit;
00212 
00213      TEUCHOS_ASSERT(Thyra::get_Epetra_Operator(*M_explicit)
00214                   ==Thyra::get_Epetra_Operator(*M_explicit_orig));
00215    
00216      //
00217      // E) Check the explicit operator
00218      //
00219    
00220      LinearOpTester<double> M_explicit_tester;
00221      M_explicit_tester.show_all_tests(true);;
00222    
00223      const bool result = M_explicit_tester.compare( *M, *M_explicit, &out );
00224      if (!result) success = false;
00225   }
00226 }
00227 
00228 } // end namespace
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines