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