Thyra Package Browser (Single Doxygen Collection) Version of the Day
test_linear_op_with_solve.cpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //    Thyra: Interfaces and Support for Abstract Numerical Algorithms
00005 //                 Copyright (2004) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include "Thyra_DefaultSerialDenseLinearOpWithSolveFactory.hpp"
00030 #include "Thyra_DefaultScaledAdjointLinearOp.hpp"
00031 #include "Thyra_DefaultAdjointLinearOpWithSolve.hpp"
00032 #include "Thyra_DefaultSpmdVectorSpace.hpp"
00033 #include "Thyra_DetachedMultiVectorView.hpp"
00034 #include "Thyra_DefaultInverseLinearOp.hpp"
00035 #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
00036 #include "Thyra_LinearOpTester.hpp"
00037 #include "Thyra_LinearOpWithSolveTester.hpp"
00038 #include "Thyra_MultiVectorStdOps.hpp"
00039 #include "Thyra_TestingTools.hpp"
00040 #include "Teuchos_GlobalMPISession.hpp"
00041 #include "Teuchos_CommandLineProcessor.hpp"
00042 #include "Teuchos_VerboseObject.hpp"
00043 #include "Teuchos_DefaultComm.hpp"
00044 #include "Teuchos_as.hpp"
00045 #include "Teuchos_StandardCatchMacros.hpp"
00046 #include "Teuchos_LocalTestingHelpers.hpp"
00047 
00048 #include "OperatorSolveHelpers.hpp"
00049 
00050 
00051 namespace Thyra {
00052 
00053 
00057 template <class Scalar>
00058 bool run_linear_op_with_solve_tests(
00059   const int n,
00060   const typename ScalarTraits<Scalar>::magnitudeType maxRelErr,
00061   const bool showAllTests,
00062   const bool dumpAll,
00063   FancyOStream &out
00064   )
00065 {
00066 
00067   using Teuchos::RCP;
00068   using Teuchos::rcp;
00069   using Teuchos::OSTab;
00070   using Teuchos::as;
00071   using Teuchos::rcp_dynamic_cast;
00072   using Teuchos::ParameterList;
00073   using Teuchos::parameterList;
00074   typedef Teuchos::ScalarTraits<Scalar> ST;
00075   typedef typename ST::magnitudeType ScalarMag;
00076   typedef Ordinal Index;
00077 
00078   out
00079     << "\n***"
00080     << "\n*** Entering run_linear_op_with_solve_tests<"<<ST::name()<<">(...) ..."
00081     << "\n***\n";
00082 
00083   bool success = true;
00084 
00085   out << "\nA) Creat a serial vector space of dimension "<<n<<" ...\n";
00086   const RCP<const Thyra::VectorSpaceBase<Scalar> >
00087     vs = Thyra::defaultSpmdVectorSpace<Scalar>(n);
00088 
00089   out << "\nB) Create a nonsingular MV object ...\n";
00090   const RCP<const MultiVectorBase<Scalar> > M =
00091     createNonsingularMultiVector(vs);
00092 
00093   out << "\nC) Create DefaultSerialDenseLinearOpWithSolve object M_lows from M ...\n";
00094 
00095   RCP<Thyra::LinearOpWithSolveBase<Scalar> >
00096     M_lows_nonconst = Thyra::linearOpWithSolve<Scalar>(
00097       *Thyra::defaultSerialDenseLinearOpWithSolveFactory<Scalar>(),
00098       M );
00099 
00100   RCP<const Thyra::LinearOpWithSolveBase<Scalar> > 
00101     M_lows = M_lows_nonconst;
00102 
00103   out << "\nD) Test the LOB interface of M_lows ...\n";
00104   Thyra::LinearOpTester<Scalar> linearOpTester;
00105   linearOpTester.set_all_error_tol(maxRelErr);
00106   linearOpTester.set_all_warning_tol(1e-2*maxRelErr);
00107   linearOpTester.show_all_tests(showAllTests);
00108   linearOpTester.dump_all(dumpAll);
00109   {
00110     OSTab tab(out);
00111     const bool result = linearOpTester.check(*M_lows, &out);
00112     if(!result) success = false;
00113   }
00114 
00115   out << "\nE) Test the LOWSB interface of M_lows ...\n";
00116   Thyra::LinearOpWithSolveTester<Scalar> linearOpWithSolveTester;
00117   {
00118     RCP<ParameterList> pl = parameterList();
00119     pl->set("All Solve Tol", maxRelErr);
00120     pl->set("All Slack Error Tol", 1e+1*maxRelErr);
00121     pl->set("All Slack Warning Tol", maxRelErr);
00122     pl->set("Show All Tests", showAllTests);
00123     pl->set("Dump All", dumpAll);
00124     linearOpWithSolveTester.setParameterList(pl);
00125   }
00126   {
00127     OSTab tab(out);
00128     const bool result = linearOpWithSolveTester.check(*M_lows, &out);
00129     if(!result) success = false;
00130   }
00131 
00132   out << "\nF) Create a DefaultInverseLinearOp object invM from M_lows and test the LOB interface ...\n";
00133   RCP<const Thyra::LinearOpBase<Scalar> > invM = inverse(M_lows);
00134   {
00135     OSTab tab(out);
00136     const bool result = linearOpTester.check(*invM, &out);
00137     if(!result) success = false;
00138   }
00139 
00140   out << "\nG) Test DefaultAdjointLinearOpWithSolve ...\n";
00141   {
00142     OSTab tab(out);
00143 
00144     out << "\nG.1) Create and test DefaultAdjointLinearOpWithSolve object M_lows_adj wrapping const M_lows ...\n";
00145     RCP<const Thyra::LinearOpWithSolveBase<Scalar> >
00146       M_lows_adj = adjointLows(M_lows);
00147 
00148     {
00149       OSTab tab2(out);
00150       
00151       out << "\nG.1.a) Test that we can extract the underlying const M_lows ...\n";
00152       {
00153         OSTab tab3(out);
00154         TEST_EQUALITY( M_lows,
00155           rcp_dynamic_cast<const Thyra::DefaultAdjointLinearOpWithSolve<Scalar> >(M_lows_adj,true)->getOp() );
00156       }
00157       
00158       out << "\nG.1.b) Testing LOB interface of DefaultAdjointLinearOpWithSolve object M_lows_adj ...\n";
00159       {
00160         OSTab tab3(out);
00161         const bool result = linearOpTester.check(*M_lows_adj, &out);
00162         if(!result) success = false;
00163       }
00164       
00165       out << "\nG.1.c) Testing LOWSB interface of DefaultAdjointLinearOpWithSolve object M_lows_adj ...\n";
00166       {
00167         OSTab tab3(out);
00168         const bool result = linearOpWithSolveTester.check(*M_lows_adj, &out);
00169         if(!result) success = false;
00170       }
00171       
00172       out << "\nG.1.d) Testing that M_lows_adj is the adjoint of M (M_adj) ...\n";
00173       const RCP<const LinearOpBase<Scalar> > M_adj = Thyra::adjoint<Scalar>(M);
00174       {
00175         OSTab tab3(out);
00176         const bool result = linearOpTester.compare(*M_lows_adj, *M_adj, &out);
00177         if(!result) success = false;
00178       }
00179 
00180     }
00181 
00182     out << "\nG.2) Create and test DefaultAdjointLinearOpWithSolve object M_lows_adj_nonconst wrapping non-const M_lows ...\n";
00183     RCP<Thyra::LinearOpWithSolveBase<Scalar> >
00184       M_lows_adj_nonconst = nonconstAdjointLows<Scalar>(M_lows_nonconst);
00185 
00186     {
00187       OSTab tab3(out);
00188       
00189       out << "\nG.2.a) Test that we can extract the underlying non-const and const M_lows ...\n";
00190       {
00191         OSTab tab4(out);
00192         TEST_EQUALITY( M_lows,
00193           rcp_dynamic_cast<Thyra::DefaultAdjointLinearOpWithSolve<Scalar> >(M_lows_adj_nonconst,true)->getOp() );
00194         TEST_EQUALITY( M_lows_nonconst,
00195           rcp_dynamic_cast<Thyra::DefaultAdjointLinearOpWithSolve<Scalar> >(M_lows_adj_nonconst,true)->getNonconstOp() );
00196       }
00197       
00198       out << "\nG.2.b) Only testing LOB interface of DefaultAdjointLinearOpWithSolve object M_lows_adj_nonconst ...\n";
00199       {
00200         OSTab tab4(out);
00201         const bool result = linearOpTester.check(*M_lows_adj_nonconst, &out);
00202         if(!result) success = false;
00203       }
00204     }
00205     
00206   }
00207   
00208   return success;
00209 
00210 }
00211 
00212 
00213 } // namespace Thyra
00214 
00215 
00216 int main( int argc, char* argv[] )
00217 {
00218 
00219   using Teuchos::CommandLineProcessor;
00220   using Teuchos::ScalarTraits;
00221   using Teuchos::as;
00222 
00223   bool success = true;
00224 
00225   Teuchos::GlobalMPISession mpiSession(&argc,&argv);
00226 
00227   Teuchos::RCP<Teuchos::FancyOStream>
00228     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00229 
00230   try {
00231 
00232     //
00233     // Read options from the command-line
00234     //
00235 
00236 
00237     CommandLineProcessor  clp;
00238     clp.throwExceptions(false);
00239     clp.addOutputSetupOptions(true);
00240 
00241     int n = 4;
00242     clp.setOption( "n", &n, "Size of the system." );
00243 
00244     double epsScale = 2e+2;
00245     clp.setOption( "eps-scale", &epsScale,
00246       "Constant (greater than 1) to scale eps by in error tests." );
00247 
00248     bool showAllTests = false;
00249     clp.setOption( "show-all-tests", "no-show-all-tests", &showAllTests,
00250       "Determines if detailed tests are shown or not." );
00251 
00252     bool dumpAll = false;
00253     clp.setOption( "dump-all", "no-dump", &dumpAll,
00254       "Determines if quantities are dumped or not." );
00255 
00256     CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
00257     if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
00258 
00259     //
00260     // Run the tests
00261     //
00262 
00263 #ifdef HAVE_THYRA_TEUCHOS_BLASFLOAT
00264     if( !Thyra::run_linear_op_with_solve_tests<float>(
00265           n, as<float>(epsScale*ScalarTraits<float>::eps()), showAllTests, dumpAll, *out)
00266       ) success = false;
00267 #endif
00268     if( !Thyra::run_linear_op_with_solve_tests<double>(
00269           n, as<double>(epsScale*ScalarTraits<double>::eps()), showAllTests, dumpAll, *out)
00270       ) success = false;
00271 #if defined(HAVE_THYRA_COMPLEX) && defined(HAVE_THYRA_TEUCHOS_BLASFLOAT)
00272     if( !Thyra::run_linear_op_with_solve_tests<std::complex<float> >(
00273           n, as<float>(epsScale*ScalarTraits<float>::eps()), showAllTests, dumpAll, *out)
00274       ) success = false;
00275 #endif
00276 #if defined(HAVE_THYRA_COMPLEX)
00277     if( !Thyra::run_linear_op_with_solve_tests<std::complex<double> >(
00278           n, as<double>(epsScale*ScalarTraits<double>::eps()), showAllTests, dumpAll, *out)
00279       ) success = false;
00280 #endif
00281 
00282   } // end try
00283   TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success);
00284 
00285   if(success)
00286     *out << "\nAll of the tests seem to have run successfully!\n";
00287   else
00288     *out << "\nOh no! at least one of the test failed!\n";  
00289   
00290   return success ? 0 : 1;
00291 
00292 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines