sillyPowerMethod_serial.cpp

Click here for a more detailed discussion of this example program.

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 "sillyPowerMethod.hpp"
00030 #include "ExampleTridiagSerialLinearOp.hpp"
00031 #include "Thyra_TestingTools.hpp"
00032 #include "Teuchos_CommandLineProcessor.hpp"
00033 #include "Teuchos_VerboseObject.hpp"
00034 #include "Teuchos_StandardCatchMacros.hpp"
00035 
00036 //
00037 // This example function is meant to show how easy it is to create
00038 // serial Thyra objects and use them with an ANA.
00039 //
00040 // This example uses a silly concrete tridiagonal matrix class
00041 // called SillyTridiagSerialLinearOp that demonstrates how
00042 // to write such subclasses.
00043 //
00044 template<class Scalar>
00045 bool runPowerMethodExample(
00046   const int    dim
00047   ,const int   maxNumIters
00048   ,const bool  verbose
00049   ,const bool  dumpAll
00050   )
00051 {
00052 
00053   using Teuchos::OSTab;
00054   typedef Teuchos::ScalarTraits<Scalar> ST;
00055   typedef typename ST::magnitudeType    ScalarMag;
00056   
00057   bool success = true;
00058   bool result;
00059 
00060   Teuchos::RCP<Teuchos::FancyOStream>
00061     out = (verbose ? Teuchos::VerboseObjectBase::getDefaultOStream() : Teuchos::null);
00062 
00063   if(verbose)
00064     *out << "\n***\n*** Running power method example using scalar type = \'" << ST::name() << "\' ...\n***\n" << std::scientific;
00065 
00066   //
00067   // (1) Setup the initial tridiagonal operator
00068   //
00069   //       [  2  -1             ]
00070   //       [ -1   2  -1         ]
00071   //  A =  [      .   .   .     ]
00072   //       [          -1  2  -1 ]
00073   //       [             -1   2 ]
00074   //
00075   if(verbose) *out << "\n(1) Constructing tridiagonal matrix A of dimension = " << dim << " ...\n";
00076   std::vector<Scalar> lower(dim-1), diag(dim), upper(dim-1);
00077   const Scalar one = ST::one(), two = Scalar(2)*one;
00078   int k = 0;
00079   diag[k] = two; upper[k] = -one;                        //  First row
00080   for( k = 1; k < dim - 1; ++k ) {
00081     lower[k-1] = -one; diag[k] = two; upper[k] = -one;   //  Middle rows
00082   }
00083   lower[k-1] = -one; diag[k] = two;                      //  Last row
00084   Teuchos::RCP<ExampleTridiagSerialLinearOp<Scalar> >
00085     A = Teuchos::rcp( new ExampleTridiagSerialLinearOp<Scalar>(dim,&lower[0],&diag[0],&upper[0]) );
00086   if( verbose && dumpAll ) *out << "\nA =\n" << *A;
00087 
00088   //
00089   // (2) Run the power method ANA
00090   //
00091   if(verbose) *out << "\n(2) Running the power method on matrix A ...\n";
00092   Scalar     lambda      = ST::zero();
00093   ScalarMag  tolerance   = ScalarMag(1e-3)*Teuchos::ScalarTraits<ScalarMag>::one();
00094   {
00095     OSTab tab(out);
00096     result = sillyPowerMethod(*A,maxNumIters,tolerance,&lambda,out.get());
00097     if(!result) success = false;
00098     if(verbose) *out << "\nEstimate of dominate eigenvalue lambda = " << lambda << std::endl;
00099   }
00100 
00101   //
00102   // (3) Increase dominance of first eigenvalue
00103   //
00104   if(verbose) *out << "\n(3) Increasing first diagonal entry by factor of 10 ...\n";
00105   diag[0] *= 10;
00106   A->initialize(dim,&lower[0],&diag[0],&upper[0]);
00107   if( verbose && dumpAll ) *out << "A =\n" << *A;
00108 
00109   //
00110   // (4) Run the power method ANA
00111   //
00112   if(verbose) *out << "\n(4) Running the power method again on matrix A ...\n";
00113   {
00114     OSTab tab(out);
00115     result = sillyPowerMethod(*A,maxNumIters,tolerance,&lambda,out.get());
00116     if(!result) success = false;
00117     if(verbose) *out << "\nEstimate of dominate eigenvalue lambda = " << lambda << std::endl;
00118   }
00119   
00120   return success;
00121 
00122 } // end runPowerMethodExample()
00123 
00124 //
00125 // Main driver program for serial implementation of the power method.
00126 //
00127 // Note that several different scalar types are used here (float,
00128 // double, std::complex<float>, std::complex<double> etc.).
00129 //
00130 int main(int argc, char *argv[])
00131 {
00132 
00133   using Teuchos::CommandLineProcessor;
00134  
00135   bool success = true;
00136   bool verbose = true;
00137   bool result;
00138 
00139   Teuchos::RCP<Teuchos::FancyOStream>
00140     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00141 
00142   try {
00143 
00144     //
00145     // Read in command-line options
00146     //
00147     
00148     int    dim          = 4;
00149     bool   dumpAll      = false;
00150 
00151     CommandLineProcessor  clp;
00152     clp.throwExceptions(false);
00153     clp.addOutputSetupOptions(true);
00154     clp.setOption( "verbose", "quiet", &verbose, "Determines if any output is printed or not." );
00155     clp.setOption( "dim", &dim, "Dimension of the linear system." );
00156     clp.setOption( "dump-all", "no-dump", &dumpAll, "Determines if quantities are dumped or not." );
00157     CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
00158     if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
00159 
00160     TEST_FOR_EXCEPTION( dim < 2, std::logic_error, "Error, dim=" << dim << " < 2 is not allowed!" );
00161 
00162     int    maxNumIters  = 10*dim;
00163     
00164     // Run using float
00165     result = runPowerMethodExample<float>(dim,maxNumIters,verbose,dumpAll);
00166     if(!result) success = false;
00167 
00168     // Run using double
00169     result = runPowerMethodExample<double>(dim,maxNumIters,verbose,dumpAll);
00170     if(!result) success = false;
00171 
00172 #if defined(HAVE_COMPLEX) && defined(HAVE_TEUCHOS_COMPLEX)
00173     
00174     // Run using std::complex<float>
00175     result = runPowerMethodExample<std::complex<float> >(dim,maxNumIters,verbose,dumpAll);
00176     if(!result) success = false;
00177     
00178     // Run using std::complex<double>
00179     result = runPowerMethodExample<std::complex<double> >(dim,maxNumIters,verbose,dumpAll);
00180     if(!result) success = false;
00181 
00182 #endif
00183 
00184 #ifdef HAVE_TEUCHOS_GNU_MP
00185     
00186     // Run using mpf_class
00187     result = runPowerMethodExample<mpf_class >(dim,maxNumIters,verbose,dumpAll);
00188     if(!result) success = false;
00189 
00190 #if defined(HAVE_COMPLEX) && defined(HAVE_TEUCHOS_COMPLEX)
00191     
00192     // Run using std::complex<mpf_class>
00193     //result = runPowerMethodExample<std::complex<mpf_class> >(dim,maxNumIters,verbose,dumpAll);
00194     //if(!result) success = false;
00195     //The above commented-out code throws a floating-point exception?
00196  
00197 #endif
00198 
00199 
00200 #endif
00201     
00202   }
00203   TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success)
00204   
00205   if (verbose) {
00206     if(success)  *out << "\nCongratulations! All of the tests checked out!\n";
00207     else         *out << "\nOh no! At least one of the tests failed!\n";
00208   }
00209   
00210   return success ? 0 : 1;
00211 
00212 } // end main()

Generated on Wed May 12 21:42:27 2010 for Thyra Operator/Vector Support by  doxygen 1.4.7