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

Generated on Thu Sep 18 12:33:02 2008 for Thyra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1