Thyra Package Browser (Single Doxygen Collection) Version of the Day
EpetraExtDiagScalingTransformer_UnitTests.cpp
Go to the documentation of this file.
00001 
00002 #include "Thyra_EpetraExtDiagScalingTransformer.hpp"
00003 #include "Thyra_DefaultScaledAdjointLinearOp.hpp"
00004 #include "Thyra_DefaultMultipliedLinearOp.hpp"
00005 #include "Thyra_DefaultDiagonalLinearOp.hpp"
00006 #include "Thyra_VectorStdOps.hpp"
00007 #include "Thyra_TestingTools.hpp"
00008 #include "Thyra_LinearOpTester.hpp"
00009 #include "Thyra_EpetraThyraWrappers.hpp"
00010 #include "Thyra_EpetraLinearOp.hpp"
00011 #include "Thyra_get_Epetra_Operator.hpp"
00012 #include "EpetraExt_readEpetraLinearSystem.h"
00013 #include "Teuchos_GlobalMPISession.hpp"
00014 #include "Teuchos_VerboseObject.hpp"
00015 #include "Teuchos_XMLParameterListHelpers.hpp"
00016 #include "Teuchos_CommandLineProcessor.hpp"
00017 #include "Teuchos_StandardCatchMacros.hpp"
00018 #ifdef _MSC_VER
00019 // This is required for operator not to be defined
00020 #include <iso646.h>
00021 #endif
00022 
00023 #ifdef HAVE_MPI
00024 #  include "Epetra_MpiComm.h"
00025 #else
00026 #  include "Epetra_SerialComm.h"
00027 #endif
00028 
00029 #include "Teuchos_UnitTestHarness.hpp"
00030 
00031 
00032 namespace {
00033 
00034 
00035 using Teuchos::null;
00036 using Teuchos::RCP;
00037 using Thyra::EpetraExtDiagScalingTransformer;
00038 using Thyra::epetraExtDiagScalingTransformer;
00039 using Thyra::VectorBase;
00040 using Thyra::LinearOpBase;
00041 using Thyra::createMember;
00042 using Thyra::LinearOpTester;
00043 using Thyra::adjoint;
00044 using Thyra::multiply;
00045 using Thyra::diagonal;
00046 
00047 
00048 std::string matrixFile = "";
00049 
00050 
00051 TEUCHOS_STATIC_SETUP()
00052 {
00053   Teuchos::UnitTestRepository::getCLP().setOption(
00054     "matrix-file", &matrixFile,
00055     "Defines the Epetra_CrsMatrix to read in."  );
00056 }
00057 
00058 // helper function to excercise all different versions of B*D*G
00059 const Teuchos::RCP<const Thyra::LinearOpBase<double> >
00060 buildADOperator(int scenario,const Teuchos::RCP<const Thyra::LinearOpBase<double> > & A,
00061                               const double vecScale, std::ostream & out)
00062 {
00063    // Create the implicit operator
00064    double scalea=-7.0;
00065    double scaled=10.0;
00066 
00067    RCP<const LinearOpBase<double> > M;
00068    RCP<VectorBase<double> > d;
00069    if(scenario<=2) 
00070       d = createMember(A->domain(), "d");
00071    else 
00072       d = createMember(A->range(), "d");
00073    V_S( d.ptr(), vecScale ); // ToDo: Set ton != 1.0 and generalize
00074  
00075    // create an operator based on requested scenario 
00076    // (these are the same as in EpetraExt)
00077    switch(scenario) {
00078    case 1: 
00079       M = multiply( scale(scalea,A), scale(scaled,diagonal(d)), "M" );
00080       out << " Testing A*D" << std::endl;
00081       break;
00082    case 2: 
00083       M = multiply( scale(scaled,diagonal(d)), scale(scalea,A), "M" );
00084       out << " Testing D*A" << std::endl;
00085       break;
00086    case 3: 
00087       M = multiply( A, scale(scaled,diagonal(d)), "M" );
00088       out << " Testing adj(A)*D" << std::endl;
00089       break;
00090    case 4: 
00091       M = multiply( scale(scaled,diagonal(d)), A, "M" );
00092       out << " Testing D*adj(A)" << std::endl;
00093       break;
00094    default:
00095       TEUCHOS_ASSERT(false);
00096       break;
00097    }
00098 
00099 
00100    out << "\nM = " << *M;
00101 
00102    return M;
00103 }
00104 
00105 TEUCHOS_UNIT_TEST( EpetraExtDiagScalingTransformer, isCompatible)
00106 {
00107    out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\'\n";
00108 
00109 #ifdef HAVE_MPI
00110    Epetra_MpiComm comm(MPI_COMM_WORLD);
00111 #else
00112    Epetra_SerialComm comm;
00113 #endif
00114 
00115    RCP<Epetra_CrsMatrix> epetra_A;
00116    RCP<Epetra_CrsMatrix> epetra_B;
00117    EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, NULL, NULL );
00118    EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_B, NULL, NULL, NULL );
00119    const RCP<const Thyra::LinearOpBase<double> > A = Thyra::epetraLinearOp(epetra_A);
00120    const RCP<const Thyra::LinearOpBase<double> > B = Thyra::epetraLinearOp(epetra_B);
00121 
00122    RCP<VectorBase<double> > d = createMember(B->domain(), "d");
00123    V_S( d.ptr(), 3.0 ); // ToDo: Set ton != 1.0 and generalize
00124 
00125    const RCP<EpetraExtDiagScalingTransformer> BD_transformer =
00126          epetraExtDiagScalingTransformer();
00127  
00128    TEST_ASSERT(BD_transformer != null);
00129 
00130    RCP<const LinearOpBase<double> > M;
00131 
00132    M = multiply(A,B);
00133    TEST_ASSERT(not BD_transformer->isCompatible(*M));
00134 
00135    M = multiply(A,scale(3.9,diagonal(d)));
00136    TEST_ASSERT(BD_transformer->isCompatible(*M));
00137 
00138    M = multiply(scale(3.0,diagonal(d)),scale(9.0,B));
00139    TEST_ASSERT(BD_transformer->isCompatible(*M));
00140 
00141    M = multiply(B,scale(3.0,diagonal(d)),scale(9.0,B));
00142    TEST_ASSERT(not BD_transformer->isCompatible(*M));
00143 }
00144 
00145 TEUCHOS_UNIT_TEST( EpetraExtDiagScalingTransformer, reapply_scaling)
00146 {
00147   //
00148   // A) Read in problem matrices
00149   //
00150   
00151   out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\'\n";
00152     
00153 #ifdef HAVE_MPI
00154   Epetra_MpiComm comm(MPI_COMM_WORLD);
00155 #else
00156   Epetra_SerialComm comm;
00157 #endif
00158   RCP<Epetra_CrsMatrix> epetra_A;
00159   EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, NULL, NULL );
00160 
00161   //
00162   // B) Create the Thyra wrapped version
00163   //
00164  
00165   const RCP<const Thyra::LinearOpBase<double> > A = Thyra::epetraLinearOp(epetra_A);
00166 
00167   RCP<const Thyra::LinearOpBase<double> > M;
00168 
00169   // test scale on right
00170   for(int scenario=1;scenario<=2;scenario++) {
00171      out << "**********************************************************" << std::endl;
00172      out << "**************** Starting Scenario = " << scenario << std::endl;
00173      out << "**********************************************************" << std::endl;
00174 
00175      //
00176      // D) Check compatibility
00177      //
00178      const RCP<EpetraExtDiagScalingTransformer> BD_transformer =
00179            epetraExtDiagScalingTransformer();
00180 
00181      TEST_ASSERT(BD_transformer != null);
00182 
00183      //
00184      // E) Do the transformation (twice)
00185      //
00186 
00187      const RCP<LinearOpBase<double> > M_explicit = BD_transformer->createOutputOp();
00188      TEST_ASSERT(M_explicit != null);
00189 
00190      // do trans
00191      M = buildADOperator(scenario,A,4.5,out);
00192      BD_transformer->transform( *M, M_explicit.ptr() );
00193      const RCP<const Epetra_Operator> eOp_explicit0 = Thyra::get_Epetra_Operator(*M_explicit);
00194 
00195      M = buildADOperator(scenario,A,7.5,out);
00196      BD_transformer->transform( *M, M_explicit.ptr() );
00197      const RCP<const Epetra_Operator> eOp_explicit1 = Thyra::get_Epetra_Operator(*M_explicit);
00198 
00199      // Epetra operator should not change!
00200      TEST_EQUALITY(eOp_explicit0,eOp_explicit1);
00201 
00202      out << "\nM_explicit = " << *M_explicit;
00203 
00204      //
00205      // F) Check the explicit operator
00206      //
00207 
00208      LinearOpTester<double> M_explicit_tester;
00209      M_explicit_tester.show_all_tests(true);;
00210 
00211      const bool result = M_explicit_tester.compare( *M, *M_explicit, &out );
00212      if (!result) success = false;
00213   }
00214 }
00215 
00216 TEUCHOS_UNIT_TEST( EpetraExtDiagScalingTransformer, basic_BDG_GDB)
00217 {
00218   
00219   //
00220   // A) Read in problem matrices
00221   //
00222   
00223   out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\'\n";
00224     
00225 #ifdef HAVE_MPI
00226   Epetra_MpiComm comm(MPI_COMM_WORLD);
00227 #else
00228   Epetra_SerialComm comm;
00229 #endif
00230   RCP<Epetra_CrsMatrix> epetra_A;
00231   EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, NULL, NULL );
00232   
00233   //
00234   // B) Create the Thyra wrapped version
00235   //
00236  
00237   const RCP<const Thyra::LinearOpBase<double> > A = Thyra::epetraLinearOp(epetra_A);
00238 
00239   //
00240   // C) Create implicit B*D*B operator
00241   //
00242 
00243   // build scenario=1 -> B*D*B, scenario=2-> B*D*B', scenario=3 -> B'*D*B
00244   for(int scenario=1;scenario<=4;scenario++) {
00245      out << "**********************************************************" << std::endl;
00246      out << "**************** Starting Scenario = " << scenario << std::endl;
00247      out << "**********************************************************" << std::endl;
00248 
00249      RCP<const Thyra::LinearOpBase<double> > M = buildADOperator(scenario,A,4.5,out);
00250 
00251      //
00252      // D) Check compatibility
00253      //
00254      const RCP<EpetraExtDiagScalingTransformer> BD_transformer =
00255            epetraExtDiagScalingTransformer();
00256 
00257      TEST_ASSERT(BD_transformer != null);
00258 
00259      BD_transformer->isCompatible(*M);
00260 
00261      //
00262      // E) Do the transformation
00263      //
00264 
00265      const RCP<LinearOpBase<double> > M_explicit = BD_transformer->createOutputOp();
00266      BD_transformer->transform( *M, M_explicit.ptr() );
00267 
00268      out << "\nM_explicit = " << *M_explicit;
00269 
00270      //
00271      // F) Check the explicit operator
00272      //
00273 
00274      LinearOpTester<double> M_explicit_tester;
00275      M_explicit_tester.show_all_tests(true);;
00276 
00277      const bool result = M_explicit_tester.compare( *M, *M_explicit, &out );
00278      if (!result) success = false;
00279   }
00280 
00281 }
00282 
00283 } // namespace
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines