Thyra Package Browser (Single Doxygen Collection) Version of the Day
test_scalar_product.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 "Teuchos_CommandLineProcessor.hpp"
00030 #include "Teuchos_VerboseObject.hpp"
00031 
00032 #ifndef SUN_CXX
00033 
00034 #include "Thyra_DefaultSpmdVectorSpace.hpp"
00035 #include "Thyra_DefaultSpmdMultiVector.hpp"
00036 #include "Thyra_LinearOpScalarProd.hpp"
00037 #include "Thyra_EuclideanScalarProd.hpp"
00038 #include "Thyra_VectorBase.hpp"
00039 #include "Thyra_MultiVectorBase.hpp"
00040 #include "Thyra_VectorStdOps.hpp"
00041 #include "Thyra_MultiVectorStdOps.hpp"
00042 #include "Thyra_DefaultDiagonalLinearOp.hpp"
00043 #include "Thyra_DefaultSerialVectorSpaceConverter.hpp"
00044 #include "Thyra_DefaultScaledAdjointLinearOp.hpp"
00045 #include "Thyra_TestingTools.hpp"
00046 #include "Thyra_LinearOpTester.hpp"
00047 
00050 template <class Scalar>
00051 bool run_scalar_product_tests(
00052   const int                                                     n
00053   ,const typename Teuchos::ScalarTraits<Scalar>::magnitudeType  &tol
00054   ,const bool                                                   dumpAll
00055   ,Teuchos::FancyOStream                                        *out_arg
00056   )
00057 {
00058 
00059   using Thyra::relErr;
00060   using Thyra::testBoolExpr;
00061   typedef Teuchos::ScalarTraits<Scalar>    ST;
00062   typedef typename ST::magnitudeType       ScalarMag;
00063   typedef Teuchos::ScalarTraits<ScalarMag> SMT;
00064   using Teuchos::Ptr;
00065   using Teuchos::RCP;
00066   using Teuchos::rcp;
00067   using Teuchos::as;
00068   using Teuchos::rcp_implicit_cast;
00069   using Teuchos::OSTab;
00070 
00071   RCP<Teuchos::FancyOStream>
00072     out = rcp(new Teuchos::FancyOStream(rcp(out_arg,false)));
00073 
00074   if(out.get()) *out << "\n*** Entering run_scalar_product_tests<"<<ST::name()<<">(...) ...\n" << std::boolalpha;
00075 
00076   bool success = true, result;
00077 
00078   RCP<Thyra::DefaultSpmdVectorSpace<Scalar> >
00079     domain = Thyra::defaultSpmdVectorSpace<Scalar>(n/2),
00080     range  = Thyra::defaultSpmdVectorSpace<Scalar>(n);
00081 
00082   RCP<Thyra::DefaultSpmdMultiVector<Scalar> >
00083     op_coeff = rcp(new Thyra::DefaultSpmdMultiVector<Scalar>(range,domain)),
00084     op       = rcp(new Thyra::DefaultSpmdMultiVector<Scalar>(range,domain));
00085 
00086   RCP<Thyra::DefaultDiagonalLinearOp<Scalar> >
00087     domainScalarProdOp = rcp(
00088       new Thyra::DefaultDiagonalLinearOp<Scalar>(
00089         rcp_implicit_cast<const Thyra::VectorSpaceBase<Scalar> >(domain)
00090         )
00091       ),
00092     rangeScalarProdOp = rcp(
00093       new Thyra::DefaultDiagonalLinearOp<Scalar>(
00094         rcp_implicit_cast<const Thyra::VectorSpaceBase<Scalar> >(range)
00095         )
00096       );
00097 
00098   Thyra::seed_randomize<Scalar>(0);
00099   Thyra::randomize<Scalar>( as<Scalar>(-ST::one()), ST::one(),
00100     Ptr<Thyra::MultiVectorBase<Scalar> >(op_coeff.ptr()) );
00101   if(out.get() && dumpAll) *out << "\nop_coeff =\n" << *op_coeff;
00102   RCP<const Thyra::VectorSpaceConverterBase<ScalarMag,Scalar> >
00103     vecSpcConverterFromMag = rcp(new Thyra::DefaultSerialVectorSpaceConverter<ScalarMag,Scalar>());
00104   RCP<const Thyra::VectorSpaceBase<ScalarMag> >
00105     magDomain = vecSpcConverterFromMag->createVectorSpaceFrom(*rcp_implicit_cast<const Thyra::VectorSpaceBase<Scalar> >(domain)),
00106     magRange  = vecSpcConverterFromMag->createVectorSpaceFrom(*rcp_implicit_cast<const Thyra::VectorSpaceBase<Scalar> >(range));
00107   RCP<Thyra::VectorBase<ScalarMag> >
00108     _domainScalarProdDiag = createMember(magDomain),
00109     _rangeScalarProdDiag  = createMember(*magRange);
00110   Thyra::randomize( ScalarMag(ScalarMag(+1)*SMT::one()), ScalarMag(ScalarMag(+2)*SMT::one()), _domainScalarProdDiag.ptr() );
00111   Thyra::randomize( ScalarMag(ScalarMag(+1)*SMT::one()), ScalarMag(ScalarMag(+2)*SMT::one()), _rangeScalarProdDiag.ptr() );
00112   vecSpcConverterFromMag->convert( *_domainScalarProdDiag, &*domainScalarProdOp->getNonconstDiag() );
00113   vecSpcConverterFromMag->convert( *_rangeScalarProdDiag, &*rangeScalarProdOp->getNonconstDiag() );
00114 
00115   RCP<const Thyra::EuclideanScalarProd<Scalar> >
00116     euclideanScalarProd = rcp(new Thyra::EuclideanScalarProd<Scalar>());
00117   RCP<const Thyra::LinearOpScalarProd<Scalar> >
00118     domainScalarProd = rcp(new Thyra::LinearOpScalarProd<Scalar>(domainScalarProdOp.create_weak())),
00119     rangeScalarProd = rcp(new Thyra::LinearOpScalarProd<Scalar>(rangeScalarProdOp.create_weak()));
00120 
00121   const ScalarMag warning_tol = ScalarMag(1e-2)*tol, error_tol = tol;
00122   Thyra::LinearOpTester<Scalar> linearOpTester;
00123   linearOpTester.linear_properties_warning_tol(warning_tol);
00124   linearOpTester.linear_properties_error_tol(error_tol);
00125   linearOpTester.adjoint_warning_tol(warning_tol);
00126   linearOpTester.adjoint_error_tol(error_tol);
00127   linearOpTester.show_all_tests(true);
00128   linearOpTester.dump_all(dumpAll);
00129 
00130   if(out.get()) *out << "\nTesting LinearOpBase with Euclidean domain and range scalar products ...\n";
00131   {
00132     OSTab tab(out);
00133     Thyra::assign( &*op, *op_coeff );
00134     if(out.get() && dumpAll) *out << "\nop =\n" << *op;
00135     if(out.get() && dumpAll) *out << "\nop' =\n" << *Thyra::adjoint(Teuchos::rcp_implicit_cast<const Thyra::LinearOpBase<Scalar> >(op));
00136     result = testBoolExpr("op->domain()->isEuclidean()",op->domain()->isEuclidean(),true,out.get());
00137     if(!result) success = false;
00138     result = testBoolExpr("op->range()->isEuclidean()",op->range()->isEuclidean(),true,out.get());
00139     if(!result) success = false;
00140     result = linearOpTester.check(*op,out.get());
00141     if(!result) success = false;
00142   }
00143  
00144   if(out.get()) *out << "\nTesting LinearOpBase with non-Euclidean domain and Euclidean range scalar products ...\n";
00145   {
00146     OSTab tab(out);
00147     range->setScalarProd(euclideanScalarProd);
00148     domain->setScalarProd(domainScalarProd);
00149     op->initialize(range,domain);
00150     Thyra::assign( &*op, *op_coeff );
00151     if(out.get() && dumpAll) *out << "\nop =\n" << *op;
00152     if(out.get() && dumpAll) *out << "\nop' =\n" << *Thyra::adjoint(Teuchos::rcp_implicit_cast<const Thyra::LinearOpBase<Scalar> >(op));
00153     result = testBoolExpr("op->domain()->isEuclidean()",op->domain()->isEuclidean(),false,out.get());
00154     if(!result) success = false;
00155     result = testBoolExpr("op->range()->isEuclidean()",op->range()->isEuclidean(),true,out.get());
00156     if(!result) success = false;
00157     result = linearOpTester.check(*op,out.get());
00158     if(!result) success = false;
00159   }
00160     
00161   if(out.get()) *out << "\nTesting LinearOpBase with Euclidean domain and non-Euclidean range scalar products ...\n";
00162   {
00163     OSTab tab(out);
00164     range->setScalarProd(rangeScalarProd);
00165     domain->setScalarProd(euclideanScalarProd);
00166     op->initialize(range,domain);
00167     Thyra::assign( &*op, *op_coeff );
00168     if(out.get() && dumpAll) *out << "\nop =\n" << *op;
00169     if(out.get() && dumpAll) *out << "\nop' =\n" << *Thyra::adjoint(Teuchos::rcp_implicit_cast<const Thyra::LinearOpBase<Scalar> >(op));
00170     result = testBoolExpr("op->domain()->isEuclidean()",op->domain()->isEuclidean(),true,out.get());
00171     if(!result) success = false;
00172     result = testBoolExpr("op->range()->isEuclidean()",op->range()->isEuclidean(),false,out.get());
00173     if(!result) success = false;
00174     result = linearOpTester.check(*op,out.get());
00175     if(!result) success = false;
00176   }
00177     
00178   if(out.get()) *out << "\nTesting LinearOpBase with non-Euclidean domain and non-Euclidean range scalar products ...\n";
00179   {
00180     OSTab tab(out);
00181     range->setScalarProd(rangeScalarProd);
00182     domain->setScalarProd(domainScalarProd);
00183     op->initialize(range,domain);
00184     Thyra::assign( &*op, *op_coeff );
00185     if(out.get() && dumpAll) *out << "\nop =\n" << *op;
00186     if(out.get() && dumpAll) *out << "\nop' =\n" << *Thyra::adjoint(Teuchos::rcp_implicit_cast<const Thyra::LinearOpBase<Scalar> >(op));
00187     result = testBoolExpr("op->domain()->isEuclidean()",op->domain()->isEuclidean(),false,out.get());
00188     if(!result) success = false;
00189     result = testBoolExpr("op->range()->isEuclidean()",op->range()->isEuclidean(),false,out.get());
00190     if(!result) success = false;
00191     result = linearOpTester.check(*op,out.get());
00192     if(!result) success = false;
00193   }
00194     
00195   if(out.get()) *out << "\n*** Leaving run_scalar_product_tests<"<<ST::name()<<">(...) ...\n";
00196 
00197   return success;
00198 
00199 } // end run_scalar_product_tests() [Doxygen looks for this!]
00200 
00201 #endif // SUN_CXX
00202 
00203 int main( int argc, char* argv[] ) {
00204 
00205   bool success = true;
00206   bool verbose = true;
00207 
00208   using Teuchos::CommandLineProcessor;
00209 
00210   Teuchos::RCP<Teuchos::FancyOStream>
00211     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00212 
00213   try {
00214 
00215     //
00216     // Read options from command-line
00217     //
00218 
00219     int n         = 4;
00220     bool dumpAll  = false;
00221 
00222     CommandLineProcessor  clp;
00223     clp.throwExceptions(false);
00224     clp.addOutputSetupOptions(true);
00225     clp.setOption( "verbose", "quiet", &verbose, "Set if output is printed or not." );
00226     clp.setOption( "n", &n, "Number of elements in each constituent vector." );
00227     clp.setOption( "dump-all", "no-dump-all", &dumpAll, "Determines if vectors are printed or not." );
00228     CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
00229     if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
00230 
00231 #ifndef SUN_CXX
00232 
00233     //
00234     // Run the tests
00235     //
00236 
00237 #ifdef HAVE_THYRA_TEUCHOS_BLASFLOAT
00238     if( !run_scalar_product_tests<float>(n,float(1e-5),dumpAll,verbose?&*out:NULL) ) success = false;
00239 #endif // HAVE_THYRA_TEUCHOS_BLASFLOAT
00240     if( !run_scalar_product_tests<double>(n,double(1e-14),dumpAll,verbose?&*out:NULL) ) success = false;
00241 #if defined(HAVE_THYRA_COMPLEX)
00242 #ifdef HAVE_THYRA_TEUCHOS_BLASFLOAT
00243     if( !run_scalar_product_tests<std::complex<float> >(n,float(1e-5),dumpAll,verbose?&*out:NULL) ) success = false;
00244 #endif // HAVE_THYRA_TEUCHOS_BLASFLOAT
00245     if( !run_scalar_product_tests<std::complex<double> >(n,double(1e-14),dumpAll,verbose?&*out:NULL) ) success = false;
00246 #endif
00247 #ifdef HAVE_TEUCHOS_GNU_MP
00248     if( !run_scalar_product_tests<mpf_class>(n,mpf_class(1e-14),dumpAll,verbose?&*out:NULL) ) success = false;
00249 #endif
00250 
00251 #endif // ifndef SUN_CXX
00252 
00253   } // end try
00254   catch( const std::exception &excpt ) {
00255     if(verbose)
00256       std::cerr << "*** Caught a standard exception : " << excpt.what() << std::endl;
00257     success = false;
00258   }
00259   catch( ... ) {
00260     if(verbose)
00261       std::cerr << "*** Caught an unknown exception!\n";
00262     success = false;
00263   }
00264 
00265 #ifndef SUN_CXX
00266 
00267   if(verbose) {
00268     if(success)
00269       *out << "\nAll of the tests seem to have run successfully!\n";
00270     else
00271       *out << "\nOh no! at least one of the test failed!\n";  
00272   }
00273   
00274   return success ? 0 : 1;
00275 
00276 #else // ifndef SUN_CXX
00277 
00278   if (verbose) {
00279     std::cout << "\nError, the test was never run since SUN_CXX was defined and this test does not build on the Sun compiler!\n";
00280   }
00281   
00282   return 1;
00283 
00284 #endif //ifndef SUN_CXX
00285 
00286 } // end main() [Doxygen looks for this!]
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines