silly1DFFT_serial.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 "ComplexFFTLinearOp.hpp"
00034 #include "RealComplexFFTLinearOp.hpp"
00035 #include "Thyra_DefaultInverseLinearOp.hpp"
00036 #include "Thyra_LinearOpTester.hpp"
00037 #include "Thyra_LinearOpWithSolveTester.hpp"
00038 #include "Thyra_ListedMultiVectorRandomizer.hpp"
00039 #include "Thyra_DefaultSerialVectorSpaceConverter.hpp"
00040 #include "Teuchos_VerboseObject.hpp"
00041 #include "Teuchos_arrayArg.hpp"
00042 #include "Teuchos_Time.hpp"
00043 #include "Teuchos_ScalarTraits.hpp"
00044 
00045 //
00046 // Creates random complex multi-vectors for FFT with symmetic entries
00047 //
00048 
00049 template<class RealScalar>
00050 class SymmetricComplexMultiVectorRandomizer : public Thyra::MultiVectorRandomizerBase< std::complex<RealScalar> > {
00051 public:
00052 
00053   typedef std::complex<RealScalar> Scalar;
00054 
00055   bool isCompatible( const Thyra::VectorSpaceBase<Scalar> &space ) const
00056     {
00057       return space.hasInCoreView();
00058     }
00059 
00060   void randomize( Thyra::MultiVectorBase<Scalar> *mv )
00061     {
00062       typedef Teuchos::ScalarTraits<Scalar> ST;
00063 #     ifdef TEUCHOS_DEBUG
00064       TEST_FOR_EXCEPT( mv == NULL );
00065       TEST_FOR_EXCEPT( mv->range()->dim() % 2 != 0 );
00066 #     endif
00067       Thyra::DetachedMultiVectorView<Scalar> ev_mv(*mv);
00068       const Thyra::Index n = ev_mv.subDim();
00069       for( Thyra::Index j = 1; j <= ev_mv.numSubCols(); ++j ) {
00070         for( Thyra::Index i = 1; i <= n/2; ++i ) {
00071           const Scalar val = ST::random();
00072           ev_mv(i,j)     = val;
00073           ev_mv(n-i+1,j) = ST::conjugate(val);
00074         }
00075       }
00076     }
00077   
00078 };
00079 
00080 //
00081 // This example program does some stuff with FFT.
00082 //
00083 template<class RealScalar>
00084 bool run1DFFTExample(
00085   const int                                                      N
00086   ,const bool                                                    verbose
00087   ,const bool                                                    dumpAll
00088   ,const RealScalar                                              tolerance
00089   ,const int                                                     outputPrec
00090   )
00091 {
00092   using Teuchos::RefCountPtr; using Teuchos::rcp;
00093   using Teuchos::OSTab;
00094   typedef std::complex<RealScalar> ComplexScalar;
00095   typedef Teuchos::ScalarTraits<RealScalar> RST;
00096   bool success = true;
00097   bool result;
00098 
00099   Teuchos::RefCountPtr<Teuchos::FancyOStream>
00100     out = ( verbose ? Teuchos::VerboseObjectBase::getDefaultOStream() : Teuchos::null );
00101 
00102   if(outputPrec > 0) out->precision(outputPrec);
00103 
00104   if(verbose)
00105     *out << "\n***\n*** Running 1D FFT example using real scalar type = \'" << RST::name() << "\' ...\n***\n";
00106 
00107   Teuchos::Time timer("");
00108   timer.start(true);
00109 
00110   if(verbose) *out << "\nA) Constructing a 1D complex-to-complex FFT linear operator C ...\n";
00111 
00112   Teuchos::RefCountPtr< const Thyra::LinearOpWithSolveBase<ComplexScalar> >
00113     C = Teuchos::rcp( new ComplexFFTLinearOp<RealScalar>(N) );
00114 
00115   if(verbose) *out << "\nB) Constructing as set of simple known vectors to be used as random domain and range vectors ...\n";
00116   Thyra::DefaultSerialVectorSpaceConverter<RealScalar,ComplexScalar>
00117     realToComplexConverter;
00118   RefCountPtr<const Thyra::VectorSpaceBase<RealScalar> >
00119     realDomainVecSpc = realToComplexConverter.createVectorSpaceFrom(*C->domain());
00120   RefCountPtr<Thyra::MultiVectorBase<RealScalar> >
00121     realDomainVec = Thyra::createMember(realDomainVecSpc);
00122   Thyra::seed_randomize<RealScalar>(0);
00123   Thyra::randomize( RealScalar(-RST::one()), RST::one(), &*realDomainVec );
00124   if(verbose && dumpAll)
00125     *out << "\nrealDomainVec:\n" << *realDomainVec;
00126   RefCountPtr<Thyra::MultiVectorBase<ComplexScalar> >
00127     complexDomainVec = Thyra::createMember(C->domain()),
00128     complexRangeVec = Thyra::createMember(C->range());
00129   realToComplexConverter.convert(*realDomainVec,&*complexDomainVec);
00130   Thyra::apply( *C, Thyra::NOTRANS, *complexDomainVec, &*complexRangeVec );
00131   if(verbose && dumpAll)
00132     *out << "\ncomplexDomainVec:\n" << *complexDomainVec << "\ncomplexRangeVec:\n" << *complexRangeVec;
00133   Thyra::ListedMultiVectorRandomizer<RealScalar>
00134     realDomainRand( Teuchos::arrayArg<RefCountPtr<const Thyra::MultiVectorBase<RealScalar> > >(realDomainVec)(), 1 );
00135   Thyra::ListedMultiVectorRandomizer<ComplexScalar>
00136     complexDomainRand( Teuchos::arrayArg<RefCountPtr<const Thyra::MultiVectorBase<ComplexScalar> > >(complexDomainVec)(), 1 ),
00137     complexRangeRand( Teuchos::arrayArg<RefCountPtr<const Thyra::MultiVectorBase<ComplexScalar> > >(complexRangeVec)(), 1 );
00138 
00139   if(verbose) *out << "\nC) Testing the LinearOpBase interface of the constructed linear operator C ...\n";
00140   Thyra::LinearOpTester<ComplexScalar> linearOpTester;
00141   linearOpTester.set_all_error_tol(tolerance);
00142   linearOpTester.set_all_warning_tol(RealScalar(RealScalar(1e-2)*tolerance));
00143   linearOpTester.show_all_tests(true);
00144   linearOpTester.dump_all(dumpAll);
00145   result = linearOpTester.check(*C,&complexRangeRand,&complexDomainRand,out.get());
00146   if(!result) success = false;
00147 
00148   if(verbose) *out << "\nD) Testing the LinearOpWithSolveBase interface of the constructed linear operator C ...\n";
00149 
00150   Thyra::LinearOpWithSolveTester<ComplexScalar> linearOpWithSolveTester;
00151   linearOpWithSolveTester.set_all_solve_tol(tolerance);
00152   linearOpWithSolveTester.set_all_slack_error_tol(RealScalar(RealScalar(1e+1)*tolerance));
00153   linearOpWithSolveTester.set_all_slack_warning_tol(tolerance);
00154   linearOpWithSolveTester.show_all_tests(true);
00155   linearOpWithSolveTester.dump_all(dumpAll);
00156   result = linearOpWithSolveTester.check(*C,out.get());
00157   if(!result) success = false;
00158 
00159   if(out.get()) *out << "\nE) Creating a DefaultInverseLinearOp object invC from C and testing the LinearOpBase interface ...\n";
00160 
00161   Teuchos::RefCountPtr<const Thyra::LinearOpBase<ComplexScalar> >
00162     invC = inverse(C);
00163 
00164   result = linearOpTester.check(*invC,out.get());
00165   if(!result) success = false;
00166 
00167   if(verbose) *out << "\nF) Constructing a 1D real-to-complex FFT linear operator R ...\n";
00168 
00169   Teuchos::RefCountPtr< const Thyra::LinearOpWithSolveBase< ComplexScalar, RealScalar > >
00170     R = Teuchos::rcp( new RealComplexFFTLinearOp<RealScalar>(N) );
00171 
00172   if(verbose) *out << "\nG) Testing the LinearOpBase interface of the constructed linear operator R ...\n";
00173   SymmetricComplexMultiVectorRandomizer<RealScalar> symmetricComplexMultiVectorRandomizer;
00174   Thyra::LinearOpTester<ComplexScalar,RealScalar> RlinearOpTester;
00175   RlinearOpTester.set_all_error_tol(tolerance);
00176   RlinearOpTester.set_all_warning_tol(RealScalar(RealScalar(1e-2)*tolerance));
00177   RlinearOpTester.show_all_tests(true);
00178   RlinearOpTester.dump_all(dumpAll);
00179   result = RlinearOpTester.check(*R,&complexRangeRand,&realDomainRand,out.get());
00180   if(!result) success = false;
00181 
00182   timer.stop();
00183 
00184   if(verbose) *out << "\nTotal time = " << timer.totalElapsedTime() << " sec\n";
00185 
00186   return success;
00187 
00188 } // end run1DFFTExample()
00189 
00190 #endif // __sun
00191 
00192 //
00193 // Actual main driver program
00194 //
00195 int main(int argc, char *argv[])
00196 {
00197 
00198   bool success = true;
00199 
00200   using Teuchos::CommandLineProcessor;
00201 
00202   bool verbose = true;
00203   bool result;
00204 
00205   Teuchos::RefCountPtr<Teuchos::FancyOStream>
00206     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00207 
00208   try {
00209 
00210     //
00211     // Read in command-line options
00212     //
00213 
00214     int    N             = 4;
00215     bool   dumpAll       = false;
00216     int    outputPrec    = -1;
00217 
00218     CommandLineProcessor  clp;
00219     clp.throwExceptions(false);
00220     clp.addOutputSetupOptions(true);
00221 
00222     clp.setOption( "verbose", "quiet", &verbose, "Determines if any output is printed or not." );
00223     clp.setOption( "N", &N, "Power of 2 for size of the FFT." );
00224     clp.setOption( "dump-all", "no-dump-all", &dumpAll, "Print all objects or not." );
00225     clp.setOption( "output-prec", &outputPrec, "Precision for outputting floating point numbers." );
00226 
00227     CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
00228     if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
00229 
00230 #ifndef __sun
00231 
00232     TEST_FOR_EXCEPTION( N < 0, std::logic_error, "Error, N=" << N << " < 1 is not allowed!" );
00233 
00234     // Run using float
00235     result = run1DFFTExample<float>(N,verbose,dumpAll,1e-5,outputPrec);
00236     if(!result) success = false;
00237 
00238     // Run using double
00239     result = run1DFFTExample<double>(N,verbose,dumpAll,1e-13,outputPrec);
00240     if(!result) success = false;
00241 
00242 #ifdef HAVE_TEUCHOS_GNU_MP
00243 
00244     // Run using mpf_class
00245     //result = run1DFFTExample<mpf_class>(N,verbose,dumpAll,1e-13,outputPrec);
00246     //if(!result) success = false;
00247 
00248 #endif // HAVE_TEUCHOS_GNU_MP
00249 
00250 #endif // ifndef __sun
00251 
00252   }
00253   catch( const std::exception &excpt ) {
00254     std::cerr << "*** Caught standard exception : " << excpt.what() << std::endl;
00255     success = false;
00256   }
00257   catch( ... ) {
00258     std::cerr << "*** Caught an unknown exception\n";
00259     success = false;
00260   }
00261 
00262 #ifndef __sun
00263 
00264   if (verbose) {
00265     if(success)   *out << "\nCongratulations! All of the tests checked out!\n";
00266     else          *out << "\nOh no! At least one of the tests failed!\n";
00267   }
00268   
00269   return success ? 0 : 1;
00270 
00271 #else // ifndef __sun
00272 
00273   if (verbose) {
00274     *out << "\nError, the test was never run since __sun was defined and this test does not build on the Sun compiler!\n";
00275   }
00276   
00277   return 1;
00278 
00279 #endif // __sun
00280 
00281 } // end main()

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