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