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 "EpetraExt_readEpetraLinearSystem.h"
00057 #include "Teuchos_GlobalMPISession.hpp"
00058 #include "Teuchos_VerboseObject.hpp"
00059 #include "Teuchos_XMLParameterListHelpers.hpp"
00060 #include "Teuchos_CommandLineProcessor.hpp"
00061 #include "Teuchos_StandardCatchMacros.hpp"
00062 
00063 #ifdef HAVE_MPI
00064 #  include "Epetra_MpiComm.h"
00065 #else
00066 #  include "Epetra_SerialComm.h"
00067 #endif
00068 
00069 #include "Teuchos_UnitTestHarness.hpp"
00070 
00071 
00072 namespace {
00073 
00074 
00075 using Teuchos::null;
00076 using Teuchos::RCP;
00077 using Thyra::EpetraExtAddTransformer;
00078 using Thyra::epetraExtAddTransformer;
00079 using Thyra::VectorBase;
00080 using Thyra::LinearOpBase;
00081 using Thyra::createMember;
00082 using Thyra::LinearOpTester;
00083 using Thyra::adjoint;
00084 using Thyra::multiply;
00085 using Thyra::diagonal;
00086 
00087 
00088 std::string matrixFile = "";
00089 std::string matrixFile2 = "";
00090 
00091 
00092 TEUCHOS_STATIC_SETUP()
00093 {
00094   Teuchos::UnitTestRepository::getCLP().setOption(
00095     "matrix-file", &matrixFile,
00096     "Defines the Epetra_CrsMatrix to read in."  );
00097   Teuchos::UnitTestRepository::getCLP().setOption(
00098     "matrix-file-2", &matrixFile2,
00099     "Defines the Epetra_CrsMatrix to read in."  );
00100 }
00101 
00102 const Teuchos::RCP<const Thyra::LinearOpBase<double> >
00103 buildAddOperator(int scenario,const Teuchos::RCP<const Thyra::LinearOpBase<double> > & A,
00104                                const Teuchos::RCP<const Thyra::LinearOpBase<double> > & B)
00105 {
00106    // build operators for the various addition/adjoint scenarios
00107    RCP<const Thyra::LinearOpBase<double> > M;
00108    
00109    switch(scenario) {
00110    case 0:
00111       M = Thyra::add(A,B,"A+B");
00112       break;
00113    case 1:
00114       M = Thyra::add(A,B,"A+adj(B)");
00115       break;
00116    case 2:
00117       M = Thyra::add(A,B,"adj(A)+B");
00118       break;
00119    case 3:
00120       M = Thyra::add(A,B,"adb(A)+adb(B)");
00121       break;
00122    default:
00123       TEUCHOS_ASSERT(false);
00124       break;
00125    }
00126 
00127    return M;
00128 }
00129 
00130 TEUCHOS_UNIT_TEST( EpetraExtAddTransformer, basic_Add )
00131 {
00132   
00133   //
00134   // A) Read in problem matrices
00135   //
00136   
00137   out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\' ...\n";
00138   out << "\nReading linear system in Epetra format from the file \'"<<matrixFile2<<"\' ...\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> epetra_A;
00146   RCP<Epetra_CrsMatrix> epetra_B;
00147   EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, NULL, NULL );
00148   EpetraExt::readEpetraLinearSystem( matrixFile2, comm, &epetra_B, NULL, NULL, NULL );
00149   
00150   //
00151   // B) Create the Thyra wrapped version
00152   //
00153   double scaleA=3.7;
00154   double scaleB=-2.9;
00155  
00156   const RCP<const Thyra::LinearOpBase<double> > A = Thyra::scale<double>(scaleA,Thyra::epetraLinearOp(epetra_B));
00157   const RCP<const Thyra::LinearOpBase<double> > B = Thyra::scale<double>(scaleB,Thyra::epetraLinearOp(epetra_B));
00158 
00159   out << "\nA = " << *A;
00160   out << "\nB = " << *B;
00161 
00162   for(int scenario=0;scenario<4;scenario++) {
00163      //
00164      // C) Create implicit A+B operator
00165      //
00166    
00167      const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,A,B);
00168    
00169      //
00170      // D) Do the transformation
00171      //
00172    
00173      const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
00174    
00175      TEST_ASSERT(ApB_transformer != null);
00176    
00177      const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
00178      ApB_transformer->transform( *M, M_explicit.ptr() );
00179    
00180      out << "\nM_explicit = " << *M_explicit;
00181    
00182      //
00183      // E) Check the explicit operator
00184      //
00185    
00186      LinearOpTester<double> M_explicit_tester;
00187      M_explicit_tester.show_all_tests(true);;
00188    
00189      const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inOutArg(out) );
00190      if (!result) success = false;
00191   }
00192 }
00193 
00194 TEUCHOS_UNIT_TEST( EpetraExtAddTransformer, mod_Add )
00195 {
00196   
00197   //
00198   // A) Read in problem matrices
00199   //
00200   
00201   out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\' ...\n";
00202   out << "\nReading linear system in Epetra format from the file \'"<<matrixFile2<<"\' ...\n";
00203     
00204 #ifdef HAVE_MPI
00205   Epetra_MpiComm comm(MPI_COMM_WORLD);
00206 #else
00207   Epetra_SerialComm comm;
00208 #endif
00209   RCP<Epetra_CrsMatrix> epetra_A;
00210   RCP<Epetra_CrsMatrix> epetra_B;
00211   EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, NULL, NULL );
00212   EpetraExt::readEpetraLinearSystem( matrixFile2, comm, &epetra_B, NULL, NULL, NULL );
00213   
00214   //
00215   // B) Create the Thyra wrapped version
00216   //
00217   double scaleA=3.7;
00218   double scaleB=-2.9;
00219  
00220   const RCP<const Thyra::LinearOpBase<double> > A = Thyra::scale<double>(scaleA,Thyra::epetraLinearOp(epetra_B));
00221   const RCP<const Thyra::LinearOpBase<double> > B = Thyra::scale<double>(scaleB,Thyra::epetraLinearOp(epetra_B));
00222 
00223   out << "\nA = " << *A;
00224   out << "\nB = " << *B;
00225 
00226   for(int scenario=0;scenario<4;scenario++) {
00227      //
00228      // C) Create implicit A+B operator
00229      //
00230    
00231      const RCP<const Thyra::LinearOpBase<double> > M = buildAddOperator(scenario,A,B);
00232    
00233      //
00234      // D) Do the transformation
00235      //
00236    
00237      const RCP<EpetraExtAddTransformer> ApB_transformer = epetraExtAddTransformer();
00238    
00239      TEST_ASSERT(ApB_transformer != null);
00240    
00241      const RCP<LinearOpBase<double> > M_explicit = ApB_transformer->createOutputOp();
00242      const RCP<LinearOpBase<double> > M_explicit_orig = M_explicit;
00243      ApB_transformer->transform( *M, M_explicit.ptr() );
00244 
00245      // do some violence to one of the operators
00246      double * view; int numEntries;
00247      epetra_B->Scale(3.2);
00248      TEUCHOS_ASSERT(epetra_B->ExtractMyRowView(3,numEntries,view)==0);
00249      for(int i=0;i<numEntries;i++) view[i] += view[i]*double(i+1);
00250 
00251      // compute multiply again
00252      ApB_transformer->transform( *M, M_explicit.ptr() );
00253 
00254      out << "\nM_explicit = " << *M_explicit;
00255 
00256      TEUCHOS_ASSERT(Thyra::get_Epetra_Operator(*M_explicit)
00257                   ==Thyra::get_Epetra_Operator(*M_explicit_orig));
00258    
00259      //
00260      // E) Check the explicit operator
00261      //
00262    
00263      LinearOpTester<double> M_explicit_tester;
00264      M_explicit_tester.show_all_tests(true);;
00265    
00266      const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inoutArg(out) );
00267      if (!result) success = false;
00268   }
00269 }
00270 
00271 } // end namespace
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines