Thyra Package Browser (Single Doxygen Collection) Version of the Day
EpetraExtDiagScalingTransformer_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_EpetraExtDiagScalingTransformer.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 "Thyra_get_Epetra_Operator.hpp"
00055 #include "EpetraExt_readEpetraLinearSystem.h"
00056 #include "Teuchos_GlobalMPISession.hpp"
00057 #include "Teuchos_VerboseObject.hpp"
00058 #include "Teuchos_XMLParameterListHelpers.hpp"
00059 #include "Teuchos_CommandLineProcessor.hpp"
00060 #include "Teuchos_StandardCatchMacros.hpp"
00061 #ifdef _MSC_VER
00062 // This is required for operator not to be defined
00063 #include <iso646.h>
00064 #endif
00065 
00066 #ifdef HAVE_MPI
00067 #  include "Epetra_MpiComm.h"
00068 #else
00069 #  include "Epetra_SerialComm.h"
00070 #endif
00071 
00072 #include "Teuchos_UnitTestHarness.hpp"
00073 
00074 
00075 namespace {
00076 
00077 
00078 using Teuchos::null;
00079 using Teuchos::RCP;
00080 using Thyra::EpetraExtDiagScalingTransformer;
00081 using Thyra::epetraExtDiagScalingTransformer;
00082 using Thyra::VectorBase;
00083 using Thyra::LinearOpBase;
00084 using Thyra::createMember;
00085 using Thyra::LinearOpTester;
00086 using Thyra::adjoint;
00087 using Thyra::multiply;
00088 using Thyra::diagonal;
00089 
00090 
00091 std::string matrixFile = "";
00092 
00093 
00094 TEUCHOS_STATIC_SETUP()
00095 {
00096   Teuchos::UnitTestRepository::getCLP().setOption(
00097     "matrix-file", &matrixFile,
00098     "Defines the Epetra_CrsMatrix to read in."  );
00099 }
00100 
00101 // helper function to excercise all different versions of B*D*G
00102 const Teuchos::RCP<const Thyra::LinearOpBase<double> >
00103 buildADOperator(int scenario,const Teuchos::RCP<const Thyra::LinearOpBase<double> > & A,
00104                               const double vecScale, std::ostream & out)
00105 {
00106    // Create the implicit operator
00107    double scalea=-7.0;
00108    double scaled=10.0;
00109 
00110    RCP<const LinearOpBase<double> > M;
00111    RCP<VectorBase<double> > d;
00112    if(scenario<=2) 
00113       d = createMember(A->domain(), "d");
00114    else 
00115       d = createMember(A->range(), "d");
00116    V_S( d.ptr(), vecScale ); // ToDo: Set ton != 1.0 and generalize
00117  
00118    // create an operator based on requested scenario 
00119    // (these are the same as in EpetraExt)
00120    switch(scenario) {
00121    case 1: 
00122       M = multiply( scale(scalea,A), scale(scaled,diagonal(d)), "M" );
00123       out << " Testing A*D" << std::endl;
00124       break;
00125    case 2: 
00126       M = multiply( scale(scaled,diagonal(d)), scale(scalea,A), "M" );
00127       out << " Testing D*A" << std::endl;
00128       break;
00129    case 3: 
00130       M = multiply( A, scale(scaled,diagonal(d)), "M" );
00131       out << " Testing adj(A)*D" << std::endl;
00132       break;
00133    case 4: 
00134       M = multiply( scale(scaled,diagonal(d)), A, "M" );
00135       out << " Testing D*adj(A)" << std::endl;
00136       break;
00137    default:
00138       TEUCHOS_ASSERT(false);
00139       break;
00140    }
00141 
00142 
00143    out << "\nM = " << *M;
00144 
00145    return M;
00146 }
00147 
00148 TEUCHOS_UNIT_TEST( EpetraExtDiagScalingTransformer, isCompatible)
00149 {
00150    out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\'\n";
00151 
00152 #ifdef HAVE_MPI
00153    Epetra_MpiComm comm(MPI_COMM_WORLD);
00154 #else
00155    Epetra_SerialComm comm;
00156 #endif
00157 
00158    RCP<Epetra_CrsMatrix> epetra_A;
00159    RCP<Epetra_CrsMatrix> epetra_B;
00160    EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, NULL, NULL );
00161    EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_B, NULL, NULL, NULL );
00162    const RCP<const Thyra::LinearOpBase<double> > A = Thyra::epetraLinearOp(epetra_A);
00163    const RCP<const Thyra::LinearOpBase<double> > B = Thyra::epetraLinearOp(epetra_B);
00164 
00165    RCP<VectorBase<double> > d = createMember(B->domain(), "d");
00166    V_S( d.ptr(), 3.0 ); // ToDo: Set ton != 1.0 and generalize
00167 
00168    const RCP<EpetraExtDiagScalingTransformer> BD_transformer =
00169          epetraExtDiagScalingTransformer();
00170  
00171    TEST_ASSERT(BD_transformer != null);
00172 
00173    RCP<const LinearOpBase<double> > M;
00174 
00175    M = multiply(A,B);
00176    TEST_ASSERT(not BD_transformer->isCompatible(*M));
00177 
00178    M = multiply(A,scale(3.9,diagonal(d)));
00179    TEST_ASSERT(BD_transformer->isCompatible(*M));
00180 
00181    M = multiply(scale(3.0,diagonal(d)),scale(9.0,B));
00182    TEST_ASSERT(BD_transformer->isCompatible(*M));
00183 
00184    M = multiply(B,scale(3.0,diagonal(d)),scale(9.0,B));
00185    TEST_ASSERT(not BD_transformer->isCompatible(*M));
00186 }
00187 
00188 TEUCHOS_UNIT_TEST( EpetraExtDiagScalingTransformer, reapply_scaling)
00189 {
00190   //
00191   // A) Read in problem matrices
00192   //
00193   
00194   out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\'\n";
00195     
00196 #ifdef HAVE_MPI
00197   Epetra_MpiComm comm(MPI_COMM_WORLD);
00198 #else
00199   Epetra_SerialComm comm;
00200 #endif
00201   RCP<Epetra_CrsMatrix> epetra_A;
00202   EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, NULL, NULL );
00203 
00204   //
00205   // B) Create the Thyra wrapped version
00206   //
00207  
00208   const RCP<const Thyra::LinearOpBase<double> > A = Thyra::epetraLinearOp(epetra_A);
00209 
00210   RCP<const Thyra::LinearOpBase<double> > M;
00211 
00212   // test scale on right
00213   for(int scenario=1;scenario<=2;scenario++) {
00214      out << "**********************************************************" << std::endl;
00215      out << "**************** Starting Scenario = " << scenario << std::endl;
00216      out << "**********************************************************" << std::endl;
00217 
00218      //
00219      // D) Check compatibility
00220      //
00221      const RCP<EpetraExtDiagScalingTransformer> BD_transformer =
00222            epetraExtDiagScalingTransformer();
00223 
00224      TEST_ASSERT(BD_transformer != null);
00225 
00226      //
00227      // E) Do the transformation (twice)
00228      //
00229 
00230      const RCP<LinearOpBase<double> > M_explicit = BD_transformer->createOutputOp();
00231      TEST_ASSERT(M_explicit != null);
00232 
00233      // do trans
00234      M = buildADOperator(scenario,A,4.5,out);
00235      BD_transformer->transform( *M, M_explicit.ptr() );
00236      const RCP<const Epetra_Operator> eOp_explicit0 = Thyra::get_Epetra_Operator(*M_explicit);
00237 
00238      M = buildADOperator(scenario,A,7.5,out);
00239      BD_transformer->transform( *M, M_explicit.ptr() );
00240      const RCP<const Epetra_Operator> eOp_explicit1 = Thyra::get_Epetra_Operator(*M_explicit);
00241 
00242      // Epetra operator should not change!
00243      TEST_EQUALITY(eOp_explicit0,eOp_explicit1);
00244 
00245      out << "\nM_explicit = " << *M_explicit;
00246 
00247      //
00248      // F) Check the explicit operator
00249      //
00250 
00251      LinearOpTester<double> M_explicit_tester;
00252      M_explicit_tester.show_all_tests(true);;
00253 
00254      const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inoutArg(out) );
00255      if (!result) success = false;
00256   }
00257 }
00258 
00259 TEUCHOS_UNIT_TEST( EpetraExtDiagScalingTransformer, basic_BDG_GDB)
00260 {
00261   
00262   //
00263   // A) Read in problem matrices
00264   //
00265   
00266   out << "\nReading linear system in Epetra format from the file \'"<<matrixFile<<"\'\n";
00267     
00268 #ifdef HAVE_MPI
00269   Epetra_MpiComm comm(MPI_COMM_WORLD);
00270 #else
00271   Epetra_SerialComm comm;
00272 #endif
00273   RCP<Epetra_CrsMatrix> epetra_A;
00274   EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A, NULL, NULL, NULL );
00275   
00276   //
00277   // B) Create the Thyra wrapped version
00278   //
00279  
00280   const RCP<const Thyra::LinearOpBase<double> > A = Thyra::epetraLinearOp(epetra_A);
00281 
00282   //
00283   // C) Create implicit B*D*B operator
00284   //
00285 
00286   // build scenario=1 -> B*D*B, scenario=2-> B*D*B', scenario=3 -> B'*D*B
00287   for(int scenario=1;scenario<=4;scenario++) {
00288      out << "**********************************************************" << std::endl;
00289      out << "**************** Starting Scenario = " << scenario << std::endl;
00290      out << "**********************************************************" << std::endl;
00291 
00292      RCP<const Thyra::LinearOpBase<double> > M = buildADOperator(scenario,A,4.5,out);
00293 
00294      //
00295      // D) Check compatibility
00296      //
00297      const RCP<EpetraExtDiagScalingTransformer> BD_transformer =
00298            epetraExtDiagScalingTransformer();
00299 
00300      TEST_ASSERT(BD_transformer != null);
00301 
00302      BD_transformer->isCompatible(*M);
00303 
00304      //
00305      // E) Do the transformation
00306      //
00307 
00308      const RCP<LinearOpBase<double> > M_explicit = BD_transformer->createOutputOp();
00309      BD_transformer->transform( *M, M_explicit.ptr() );
00310 
00311      out << "\nM_explicit = " << *M_explicit;
00312 
00313      //
00314      // F) Check the explicit operator
00315      //
00316 
00317      LinearOpTester<double> M_explicit_tester;
00318      M_explicit_tester.show_all_tests(true);
00319 
00320      const bool result = M_explicit_tester.compare( *M, *M_explicit, Teuchos::inoutArg(out) );
00321      if (!result) success = false;
00322   }
00323 
00324 }
00325 
00326 } // namespace
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines