Thyra Package Browser (Single Doxygen Collection) Version of the Day
EpetraExtDiagScaledMatProdTransformer_UnitTests.cpp
Go to the documentation of this file.
00001 
00002 #include "Thyra_EpetraExtDiagScaledMatProdTransformer.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 "EpetraExt_readEpetraLinearSystem.h"
00012 #include "Teuchos_GlobalMPISession.hpp"
00013 #include "Teuchos_VerboseObject.hpp"
00014 #include "Teuchos_XMLParameterListHelpers.hpp"
00015 #include "Teuchos_CommandLineProcessor.hpp"
00016 #include "Teuchos_StandardCatchMacros.hpp"
00017 
00018 #ifdef HAVE_MPI
00019 #  include "Epetra_MpiComm.h"
00020 #else
00021 #  include "Epetra_SerialComm.h"
00022 #endif
00023 
00024 #include "Teuchos_UnitTestHarness.hpp"
00025 
00026 
00027 namespace {
00028 
00029 
00030 using Teuchos::null;
00031 using Teuchos::RCP;
00032 using Thyra::EpetraExtDiagScaledMatProdTransformer;
00033 using Thyra::epetraExtDiagScaledMatProdTransformer;
00034 using Thyra::VectorBase;
00035 using Thyra::LinearOpBase;
00036 using Thyra::createMember;
00037 using Thyra::LinearOpTester;
00038 using Thyra::adjoint;
00039 using Thyra::multiply;
00040 using Thyra::diagonal;
00041 
00042 
00043 std::string matrixFile = "";
00044 std::string matrixFile2 = "";
00045 
00046 
00047 TEUCHOS_STATIC_SETUP()
00048 {
00049   Teuchos::UnitTestRepository::getCLP().setOption(
00050     "matrix-file", &matrixFile,
00051     "Defines the Epetra_CrsMatrix to read in."  );
00052   Teuchos::UnitTestRepository::getCLP().setOption(
00053     "matrix-file-2", &matrixFile2,
00054     "Defines the second Epetra_CrsMatrix to read in."  );
00055 }
00056 
00057 // helper function to excercise all different versions of B*D*G
00058 const Teuchos::RCP<const Thyra::LinearOpBase<double> >
00059 buildBDGOperator(int scenario,const Teuchos::RCP<const Thyra::LinearOpBase<double> > & B,
00060                               const Teuchos::RCP<const Thyra::LinearOpBase<double> > & G,
00061                               const double vecScale, std::ostream & out)
00062 {
00063    // Create the implicit operator
00064    double scalea=10.0;
00065    double scaleb=-7.0;
00066    double scaled=52.0;
00067 
00068    RCP<const LinearOpBase<double> > M;
00069    RCP<VectorBase<double> > d;
00070    if(scenario<=2) 
00071       d = createMember(B->domain(), "d");
00072    else 
00073       d = createMember(B->range(), "d");
00074    V_S( d.ptr(), vecScale ); // ToDo: Set ton != 1.0 and generalize
00075  
00076    // create an operator based on requested scenario 
00077    // (these are the same as in EpetraExt)
00078    switch(scenario) {
00079    case 1: 
00080       M = multiply( scale(scalea,B), diagonal(d), scale(scaleb,G), "M" );
00081       out << " Testing B*D*G" << std::endl;
00082       break;
00083    case 2: 
00084       M = multiply( scale(scalea,B), diagonal(d), adjoint(G), "M" );
00085       out << " Testing B*D*adj(G)" << std::endl;
00086       break;
00087    case 3: 
00088       M = multiply( adjoint(B), diagonal(d), scale(scaleb,G), "M" );
00089       out << " Testing adj(B)*D*G" << std::endl;
00090       break;
00091    case 4: 
00092       M = multiply( adjoint(B), diagonal(d), adjoint(G), "M" );
00093       out << " Testing adj(B)*D*adj(G)" << std::endl;
00094       break;
00095    case 5: 
00096       M = multiply( B, scale(scaled,diagonal(d)), G, "M" );
00097       out << " Testing B*(52.0*D)*G" << std::endl;
00098       break;
00099    default:
00100       TEUCHOS_ASSERT(false);
00101       break;
00102    }
00103 
00104 
00105    out << "\nM = " << *M;
00106 
00107    return M;
00108 }
00109 
00110 // helper function to excercise all different versions of B*D*G
00111 const Teuchos::RCP<const Thyra::LinearOpBase<double> >
00112 buildBGOperator(int scenario,const Teuchos::RCP<const Thyra::LinearOpBase<double> > & B,
00113                              const Teuchos::RCP<const Thyra::LinearOpBase<double> > & G,
00114                              std::ostream & out)
00115 {
00116    // Create the implicit operator
00117    double scalea=10.0;
00118    double scaleb=-7.0;
00119    RCP<const LinearOpBase<double> > M;
00120 
00121    // create an operator based on requested scenario 
00122    // (these are the same as in EpetraExt)
00123    switch(scenario) {
00124    case 1: 
00125       M = multiply( scale(scalea,B), scale(scaleb,G), "M" );
00126       out << " Testing B*G" << std::endl;
00127       break;
00128    case 2: 
00129       M = multiply( scale(scalea,B), adjoint(G), "M" );
00130       out << " Testing B*adj(G)" << std::endl;
00131       break;
00132    case 3: 
00133       M = multiply( adjoint(B), scale(scaleb,G), "M" );
00134       out << " Testing adj(B)*G" << std::endl;
00135       break;
00136    case 4: 
00137       M = multiply( adjoint(B), adjoint(G), "M" );
00138       out << " Testing adj(B)*adj(G)" << std::endl;
00139       break;
00140    default:
00141       TEUCHOS_ASSERT(false);
00142       break;
00143    }
00144 
00145 
00146    out << "\nM = " << *M;
00147 
00148    return M;
00149 }
00150 
00151 const Teuchos::RCP<const Thyra::LinearOpBase<double> >
00152 buildBDBOperator(int scenario,const Teuchos::RCP<const Thyra::LinearOpBase<double> > & B,
00153                  const double vecScale, std::ostream & out)
00154 {
00155    return buildBDGOperator(scenario,B,B,vecScale,out);
00156 }
00157 
00158 
00159 
00160 TEUCHOS_UNIT_TEST( EpetraExtDiagScaledMatProdTransformer, basic_BDB )
00161 {
00162   
00163   //
00164   // A) Read in problem matrices
00165   //
00166   
00167   out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\' ...\n";
00168     
00169 #ifdef HAVE_MPI
00170   Epetra_MpiComm comm(MPI_COMM_WORLD);
00171 #else
00172   Epetra_SerialComm comm;
00173 #endif
00174   RCP<Epetra_CrsMatrix> epetra_B;
00175   EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_B, NULL, NULL, NULL );
00176   
00177   //
00178   // B) Create the Thyra wrapped version
00179   //
00180  
00181   const RCP<const Thyra::LinearOpBase<double> > B = Thyra::epetraLinearOp(epetra_B);
00182 
00183   //
00184   // C) Create implicit B*D*B operator
00185   //
00186 
00187   // build scenario=1 -> B*D*B, scenario=2 -> B*D*B',
00188   //       scenario=3 -> B'*D*B, scenario=4 -> B'*D*B'
00189   //int scenario = 3;
00190   for(int scenario=1;scenario<=5;scenario++) {
00191      const RCP<const Thyra::LinearOpBase<double> > M = buildBDBOperator(scenario,B,4.5,out);
00192 
00193      //
00194      // D) Do the transformation
00195      //
00196 
00197      const RCP<EpetraExtDiagScaledMatProdTransformer> BtDB_transformer =
00198        epetraExtDiagScaledMatProdTransformer();
00199 
00200      TEST_ASSERT(BtDB_transformer != null);
00201 
00202      const RCP<LinearOpBase<double> > M_explicit = BtDB_transformer->createOutputOp();
00203 
00204      BtDB_transformer->transform( *M, M_explicit.ptr() );
00205 
00206      out << "\nM_explicit = " << *M_explicit;
00207 
00208      //
00209      // E) Check the explicit operator
00210      //
00211 
00212      LinearOpTester<double> M_explicit_tester;
00213      M_explicit_tester.show_all_tests(true);;
00214 
00215      const bool result = M_explicit_tester.compare( *M, *M_explicit, &out );
00216      if (!result) success = false;
00217   }
00218 
00219 }
00220 
00221 TEUCHOS_UNIT_TEST( EpetraExtDiagScaledMatProdTransformer, basic_BDG_GDB)
00222 {
00223   
00224   //
00225   // A) Read in problem matrices
00226   //
00227   
00228   out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\'";
00229   out << " and from the file \'"<<matrixFile2<<"\' ...\n";
00230     
00231 #ifdef HAVE_MPI
00232   Epetra_MpiComm comm(MPI_COMM_WORLD);
00233 #else
00234   Epetra_SerialComm comm;
00235 #endif
00236   RCP<Epetra_CrsMatrix> epetra_B;
00237   RCP<Epetra_CrsMatrix> epetra_G;
00238   EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_B, NULL, NULL, NULL );
00239   EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_G, NULL, NULL, NULL );
00240   
00241   //
00242   // B) Create the Thyra wrapped version
00243   //
00244  
00245   const RCP<const Thyra::LinearOpBase<double> > B = Thyra::epetraLinearOp(epetra_B);
00246   const RCP<const Thyra::LinearOpBase<double> > G = Thyra::epetraLinearOp(epetra_G);
00247 
00248   //
00249   // C) Create implicit B*D*B operator
00250   //
00251 
00252   // build scenario=1 -> B*D*B, scenario=2-> B*D*B', scenario=3 -> B'*D*B
00253   for(int scenario=1;scenario<=3;scenario++) {
00254      RCP<const Thyra::LinearOpBase<double> > M;
00255      if(scenario==1 || scenario==3)
00256         M = buildBDGOperator(scenario,B,G,4.5,out);
00257      else
00258         M = buildBDGOperator(scenario,G,B,4.5,out);
00259 
00260      //
00261      // D) Do the transformation
00262      //
00263 
00264      const RCP<EpetraExtDiagScaledMatProdTransformer> BtDB_transformer =
00265        epetraExtDiagScaledMatProdTransformer();
00266 
00267      TEST_ASSERT(BtDB_transformer != null);
00268 
00269      const RCP<LinearOpBase<double> > M_explicit = BtDB_transformer->createOutputOp();
00270 
00271      BtDB_transformer->transform( *M, M_explicit.ptr() );
00272 
00273      out << "\nM_explicit = " << *M_explicit;
00274 
00275      //
00276      // E) Check the explicit operator
00277      //
00278 
00279      LinearOpTester<double> M_explicit_tester;
00280      M_explicit_tester.show_all_tests(true);;
00281 
00282      const bool result = M_explicit_tester.compare( *M, *M_explicit, &out );
00283      if (!result) success = false;
00284   }
00285 
00286 }
00287 
00288 TEUCHOS_UNIT_TEST( EpetraExtDiagScaledMatProdTransformer, basic_GDG )
00289 {
00290   
00291   //
00292   // A) Read in problem matrices
00293   //
00294   
00295   out << "\nReading linear system in Epetra format from the file \'"<<matrixFile2<<"\' ...\n";
00296     
00297 #ifdef HAVE_MPI
00298   Epetra_MpiComm comm(MPI_COMM_WORLD);
00299 #else
00300   Epetra_SerialComm comm;
00301 #endif
00302   RCP<Epetra_CrsMatrix> epetra_G;
00303   EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_G, NULL, NULL, NULL );
00304   
00305   //
00306   // B) Create the Thyra wrapped version
00307   //
00308  
00309   const RCP<const Thyra::LinearOpBase<double> > G = Thyra::epetraLinearOp(epetra_G);
00310 
00311   //
00312   // C) Create implicit B*D*B operator
00313   //
00314 
00315   // build scenario=1 -> B*D*B, scenario=3 -> B'*D*B
00316   int scenes[] = {1,4};
00317   for(int i=0;i<2;i++) {
00318      int scenario = scenes[i];
00319      const RCP<const Thyra::LinearOpBase<double> > M = buildBDGOperator(scenario,G,G,4.5,out);
00320 
00321      //
00322      // D) Do the transformation
00323      //
00324 
00325      const RCP<EpetraExtDiagScaledMatProdTransformer> BtDB_transformer =
00326        epetraExtDiagScaledMatProdTransformer();
00327 
00328      TEST_ASSERT(BtDB_transformer != null);
00329 
00330      const RCP<LinearOpBase<double> > M_explicit = BtDB_transformer->createOutputOp();
00331 
00332      BtDB_transformer->transform( *M, M_explicit.ptr() );
00333 
00334      out << "\nM_explicit = " << *M_explicit;
00335 
00336      //
00337      // E) Check the explicit operator
00338      //
00339 
00340      LinearOpTester<double> M_explicit_tester;
00341      M_explicit_tester.show_all_tests(true);;
00342 
00343      const bool result = M_explicit_tester.compare( *M, *M_explicit, &out );
00344      if (!result) success = false;
00345   }
00346 
00347 }
00348 
00349 TEUCHOS_UNIT_TEST( EpetraExtDiagScaledMatProdTransformer, basic_BG_GB_GG)
00350 {
00351   
00352   //
00353   // A) Read in problem matrices
00354   //
00355   
00356   out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\'";
00357   out << " and from the file \'"<<matrixFile2<<"\' ...\n";
00358     
00359 #ifdef HAVE_MPI
00360   Epetra_MpiComm comm(MPI_COMM_WORLD);
00361 #else
00362   Epetra_SerialComm comm;
00363 #endif
00364   RCP<Epetra_CrsMatrix> epetra_B;
00365   RCP<Epetra_CrsMatrix> epetra_G;
00366   EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_B, NULL, NULL, NULL );
00367   EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_G, NULL, NULL, NULL );
00368   
00369   //
00370   // B) Create the Thyra wrapped version
00371   //
00372  
00373   const RCP<const Thyra::LinearOpBase<double> > B = Thyra::epetraLinearOp(epetra_B);
00374   const RCP<const Thyra::LinearOpBase<double> > G = Thyra::epetraLinearOp(epetra_G);
00375 
00376   //
00377   // C) Create implicit B*D*B operator
00378   //
00379 
00380   // build scenario=1 -> B*D*B, scenario=2-> B*D*B', scenario=3 -> B'*D*B
00381   for(int scenario=1;scenario<=4;scenario++) {
00382      RCP<const Thyra::LinearOpBase<double> > M;
00383      if(scenario==1 || scenario==3)
00384         M = buildBGOperator(scenario,B,G,out);
00385      else if(scenario==2)
00386         M = buildBGOperator(scenario,G,B,out);
00387      else
00388         M = buildBGOperator(scenario,G,G,out);
00389 
00390      //
00391      // D) Do the transformation
00392      //
00393 
00394      const RCP<EpetraExtDiagScaledMatProdTransformer> BtDB_transformer =
00395        epetraExtDiagScaledMatProdTransformer();
00396 
00397      TEST_ASSERT(BtDB_transformer != null);
00398 
00399      const RCP<LinearOpBase<double> > M_explicit = BtDB_transformer->createOutputOp();
00400 
00401      BtDB_transformer->transform( *M, M_explicit.ptr() );
00402 
00403      out << "\nM_explicit = " << *M_explicit;
00404 
00405      //
00406      // E) Check the explicit operator
00407      //
00408 
00409      LinearOpTester<double> M_explicit_tester;
00410      M_explicit_tester.show_all_tests(true);;
00411 
00412      const bool result = M_explicit_tester.compare( *M, *M_explicit, &out );
00413      if (!result) success = false;
00414   }
00415 
00416 }
00417 
00418 } // namespace
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines