Thyra Package Browser (Single Doxygen Collection) Version of the Day
sillyPowerMethod_epetra.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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #include "sillyPowerMethod.hpp"
00043 #include "createTridiagEpetraLinearOp.hpp"
00044 #include "Thyra_TestingTools.hpp"
00045 #include "Thyra_EpetraLinearOp.hpp"
00046 #include "Thyra_get_Epetra_Operator.hpp"
00047 #include "Teuchos_GlobalMPISession.hpp"
00048 #include "Teuchos_CommandLineProcessor.hpp"
00049 #include "Teuchos_StandardCatchMacros.hpp"
00050 #include "Teuchos_VerboseObject.hpp"
00051 #include "Teuchos_dyn_cast.hpp"
00052 #include "Epetra_CrsMatrix.h"
00053 #include "Epetra_Map.h"
00054 
00055 //
00056 // Increase the first diagonal element of your tridiagonal matrix.
00057 //
00058 void scaleFirstDiagElement( const double diagScale, Thyra::LinearOpBase<double> *A )
00059 {
00060   using Teuchos::RCP;
00061   TEUCHOS_TEST_FOR_EXCEPT(A==NULL);
00062   // (A) Get at the underlying Epetra_Operator object that the EpetraLinearOp
00063   // object directly maintains.
00064   const RCP<Epetra_Operator> epetra_op = Thyra::get_Epetra_Operator(*A);
00065   // (B) Perform a dynamic cast to Epetra_CrsMatrix.
00066   // Note, the dyn_cast<>() template function will throw std::bad_cast
00067   // with a nice error message if the cast fails! 
00068   Epetra_CrsMatrix &crsMatrix = Teuchos::dyn_cast<Epetra_CrsMatrix>(*epetra_op);
00069   // (C) Query if the first row of the matrix is on this processor and if
00070   // it is get it and scale it.
00071   if(crsMatrix.MyGlobalRow(0)) {
00072     const int numRowNz = crsMatrix.NumGlobalEntries(0);
00073     TEUCHOS_TEST_FOR_EXCEPT( numRowNz != 2 );
00074     int returndNumRowNz; double rowValues[2]; int rowIndexes[2];
00075     crsMatrix.ExtractGlobalRowCopy( 0, numRowNz, returndNumRowNz, &rowValues[0], &rowIndexes[0] );
00076     TEUCHOS_TEST_FOR_EXCEPT( returndNumRowNz != 2 );
00077     for( int k = 0; k < numRowNz; ++k) if (rowIndexes[k] == 0) rowValues[k] *= diagScale;
00078     crsMatrix.ReplaceGlobalValues( 0, numRowNz, rowValues, rowIndexes );
00079   }
00080 } // end scaleFirstDiagElement()
00081 
00082 //
00083 // Main driver program for epetra implementation of the power method.
00084 //
00085 int main(int argc, char *argv[])
00086 {
00087 
00088   using Teuchos::outArg;
00089   using Teuchos::CommandLineProcessor;
00090   using Teuchos::RCP;
00091  
00092   bool success = true;
00093   bool result;
00094   
00095   //
00096   // (A) Setup and get basic MPI info
00097   //
00098   
00099   Teuchos::GlobalMPISession mpiSession(&argc,&argv);
00100 #ifdef HAVE_MPI
00101   MPI_Comm mpiComm = MPI_COMM_WORLD;
00102 #endif
00103   
00104   //
00105   // (B) Setup the output stream (do output only on root process!)
00106   //
00107 
00108   Teuchos::RCP<Teuchos::FancyOStream>
00109     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00110   
00111   try {
00112 
00113     //
00114     // (C) Read in commandline options
00115     //
00116 
00117 
00118     CommandLineProcessor  clp;
00119     clp.throwExceptions(false);
00120     clp.addOutputSetupOptions(true);
00121 
00122     int globalDim = 500;
00123     clp.setOption( "global-dim", &globalDim, "Global dimension of the linear system." );
00124 
00125     bool dumpAll = false;
00126     clp.setOption( "dump-all", "no-dump", &dumpAll, "Determines if quantities are dumped or not." );
00127 
00128     CommandLineProcessor::EParseCommandLineReturn
00129       parse_return = clp.parse(argc,argv);
00130     if (parse_return != CommandLineProcessor::PARSE_SUCCESSFUL)
00131       return parse_return;
00132 
00133     TEUCHOS_TEST_FOR_EXCEPTION( globalDim < 2, std::logic_error, "Error, globalDim=" << globalDim << " < 2 is not allowed!" );
00134 
00135     *out << "\n***\n*** Running power method example using Epetra implementation\n***\n" << std::scientific;
00136 
00137     //
00138     // (D) Setup the operator and run the power method!
00139     //
00140 
00141     //
00142     // (1) Setup the initial tridiagonal operator
00143     //
00144     //       [  2  -1             ]
00145     //       [ -1   2  -1         ]
00146     //  A =  [      .   .   .     ]
00147     //       [          -1  2  -1 ]
00148     //       [             -1   2 ]
00149     //
00150     *out << "\n(1) Constructing tridagonal Epetra matrix A of global dimension = " << globalDim << " ...\n";
00151     RCP<Epetra_Operator>
00152       A_epetra = createTridiagEpetraLinearOp(
00153         globalDim,
00154 #ifdef HAVE_MPI
00155         mpiComm,
00156 #endif
00157         1.0, true, *out
00158         );
00159     // Wrap in an Thyra::EpetraLinearOp object
00160     RCP<Thyra::LinearOpBase<double> >
00161       A = Thyra::nonconstEpetraLinearOp(A_epetra);
00162     //
00163     if (dumpAll) *out << "\nA =\n" << *A; // This works even in parallel!
00164     
00165     //
00166     // (2) Run the power method ANA
00167     //
00168     *out << "\n(2) Running the power method on matrix A ...\n";
00169     double  lambda      = 0.0;
00170     double  tolerance   = 1e-3;
00171     int     maxNumIters = 10*globalDim;
00172     result = sillyPowerMethod<double>(*A, maxNumIters, tolerance, outArg(lambda), *out);
00173     if(!result) success = false;
00174     *out << "\n  Estimate of dominate eigenvalue lambda = " << lambda << std::endl;
00175     
00176     //
00177     // (3) Increase dominance of first eigenvalue
00178     //
00179     *out << "\n(3) Scale the diagonal of A by a factor of 10 ...\n";
00180     scaleFirstDiagElement( 10.0, &*A );
00181     
00182     //
00183     // (4) Run the power method ANA again
00184     //
00185     *out << "\n(4) Running the power method again on matrix A ...\n";
00186     result = sillyPowerMethod<double>(*A, maxNumIters, tolerance, outArg(lambda), *out);
00187     if(!result) success = false;
00188     *out << "\n  Estimate of dominate eigenvalue lambda = " << lambda << std::endl;
00189     
00190   }
00191   TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success)
00192   
00193   if (success)
00194     *out << "\nCongratulations! All of the tests checked out!\n";
00195   else
00196     *out << "\nOh no! At least one of the tests failed!\n";
00197 
00198   return success ? 0 : 1;
00199 
00200 } // end main()
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines