sillyPowerMethod_epetra.cpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //                Epetra: Linear Algebra Services Package 
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 "createTridiagEpetraLinearOp.hpp"
00031 #include "Thyra_TestingTools.hpp"
00032 #include "Thyra_EpetraLinearOp.hpp"
00033 #include "Thyra_get_Epetra_Operator.hpp"
00034 #include "Teuchos_GlobalMPISession.hpp"
00035 #include "Teuchos_CommandLineProcessor.hpp"
00036 #include "Teuchos_StandardCatchMacros.hpp"
00037 #include "Teuchos_VerboseObject.hpp"
00038 #include "Teuchos_dyn_cast.hpp"
00039 #include "Epetra_CrsMatrix.h"
00040 #include "Epetra_Map.h"
00041 
00042 //
00043 // Increase the first diagonal element of your tridiagonal matrix.
00044 //
00045 void scaleFirstDiagElement( const double diagScale, Thyra::LinearOpBase<double> *A )
00046 {
00047   using Teuchos::RefCountPtr;
00048   TEST_FOR_EXCEPT(A==NULL);
00049   // (A) Get at the underlying Epetra_Operator object that the EpetraLinearOp
00050   // object directly maintains.
00051   const RefCountPtr<Epetra_Operator> epetra_op = Thyra::get_Epetra_Operator(*A);
00052   // (B) Perform a dynamic cast to Epetra_CrsMatrix.
00053   // Note, the dyn_cast<>() template function will throw std::bad_cast
00054   // with a nice error message if the cast fails! 
00055   Epetra_CrsMatrix &crsMatrix = Teuchos::dyn_cast<Epetra_CrsMatrix>(*epetra_op);
00056   // (C) Query if the first row of the matrix is on this processor and if
00057   // it is get it and scale it.
00058   if(crsMatrix.MyGlobalRow(0)) {
00059     const int numRowNz = crsMatrix.NumGlobalEntries(0);
00060     TEST_FOR_EXCEPT( numRowNz != 2 );
00061     int returndNumRowNz; double rowValues[2]; int rowIndexes[2];
00062     crsMatrix.ExtractGlobalRowCopy( 0, numRowNz, returndNumRowNz, &rowValues[0], &rowIndexes[0] );
00063     TEST_FOR_EXCEPT( returndNumRowNz != 2 );
00064     for( int k = 0; k < numRowNz; ++k) if (rowIndexes[k] == 0) rowValues[k] *= diagScale;
00065     crsMatrix.ReplaceGlobalValues( 0, numRowNz, rowValues, rowIndexes );
00066   }
00067 } // end scaleFirstDiagElement()
00068 
00069 //
00070 // Main driver program for epetra implementation of the power method.
00071 //
00072 int main(int argc, char *argv[])
00073 {
00074   using Teuchos::CommandLineProcessor;
00075   using Teuchos::RefCountPtr;
00076  
00077   bool success = true;
00078   bool verbose = true;
00079   bool result;
00080   
00081   //
00082   // (A) Setup and get basic MPI info
00083   //
00084   
00085   Teuchos::GlobalMPISession mpiSession(&argc,&argv);
00086   const int procRank = Teuchos::GlobalMPISession::getRank();
00087   const int numProc  = Teuchos::GlobalMPISession::getNProc();
00088 #ifdef HAVE_MPI
00089   MPI_Comm mpiComm = MPI_COMM_WORLD;
00090 #endif
00091   
00092   //
00093   // (B) Setup the output stream (do output only on root process!)
00094   //
00095 
00096   Teuchos::RefCountPtr<Teuchos::FancyOStream>
00097     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00098   
00099   try {
00100 
00101     //
00102     // (C) Read in commandline options
00103     //
00104 
00105     int    globalDim                  = 500;
00106     bool   dumpAll                    = false;
00107 
00108     CommandLineProcessor  clp;
00109     clp.throwExceptions(false);
00110     clp.addOutputSetupOptions(true);
00111 
00112     clp.setOption( "verbose", "quiet", &verbose, "Determines if any output is printed or not." );
00113     clp.setOption( "global-dim", &globalDim, "Global dimension of the linear system." );
00114     clp.setOption( "dump-all", "no-dump", &dumpAll, "Determines if quantities are dumped or not." );
00115 
00116     CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
00117     if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
00118 
00119     TEST_FOR_EXCEPTION( globalDim < 2, std::logic_error, "Error, globalDim=" << globalDim << " < 2 is not allowed!" );
00120 
00121     if(verbose) *out << "\n***\n*** Running power method example using Epetra implementation\n***\n" << std::scientific;
00122 
00123     //
00124     // (D) Setup the operator and run the power method!
00125     //
00126 
00127     //
00128     // (1) Setup the initial tridiagonal operator
00129     //
00130     //       [  2  -1             ]
00131     //       [ -1   2  -1         ]
00132     //  A =  [      .   .   .     ]
00133     //       [          -1  2  -1 ]
00134     //       [             -1   2 ]
00135     //
00136     if(verbose) *out << "\n(1) Constructing tridagonal Epetra matrix A of global dimension = " << globalDim << " ...\n";
00137     RefCountPtr<Epetra_Operator>
00138       A_epetra = createTridiagEpetraLinearOp(
00139         globalDim
00140 #ifdef HAVE_MPI
00141         ,mpiComm
00142 #endif
00143         ,1.0,verbose,*out
00144         );
00145     // Wrap in an Thyra::EpetraLinearOp object
00146     RefCountPtr<Thyra::LinearOpBase<double> >
00147       A = rcp(new Thyra::EpetraLinearOp(A_epetra));
00148     //
00149     if( verbose && dumpAll ) *out << "\nA =\n" << *A; // This works even in parallel!
00150     
00151     //
00152     // (2) Run the power method ANA
00153     //
00154     if(verbose) *out << "\n(2) Running the power method on matrix A ...\n";
00155     double  lambda      = 0.0;
00156     double  tolerance   = 1e-3;
00157     int     maxNumIters = 10*globalDim;
00158     result = sillyPowerMethod<double>(*A,maxNumIters,tolerance,&lambda,(verbose?&*out:NULL));
00159     if(!result) success = false;
00160     if(verbose) *out << "\n  Estimate of dominate eigenvalue lambda = " << lambda << std::endl;
00161     
00162     //
00163     // (3) Increase dominance of first eigenvalue
00164     //
00165     if(verbose) *out << "\n(3) Scale the diagonal of A by a factor of 10 ...\n";
00166     scaleFirstDiagElement( 10.0, &*A );
00167     
00168     //
00169     // (4) Run the power method ANA again
00170     //
00171     if(verbose) *out << "\n(4) Running the power method again on matrix A ...\n";
00172     result = sillyPowerMethod<double>(*A,maxNumIters,tolerance,&lambda,(verbose?&*out:NULL));
00173     if(!result) success = false;
00174     if(verbose) *out << "\n  Estimate of dominate eigenvalue lambda = " << lambda << std::endl;
00175     
00176   }
00177   TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success)
00178   
00179   if (verbose) {
00180     if(success)  *out << "\nCongratulations! All of the tests checked out!\n";
00181     else         *out << "\nOh no! At least one of the tests failed!\n";
00182   }
00183 
00184   return success ? 0 : 1;
00185 
00186 } // end main()

Generated on Thu Sep 18 12:37:58 2008 for Epetra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1