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 //                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::RCP;
00048   TEST_FOR_EXCEPT(A==NULL);
00049   // (A) Get at the underlying Epetra_Operator object that the EpetraLinearOp
00050   // object directly maintains.
00051   const RCP<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 
00075   using Teuchos::outArg;
00076   using Teuchos::CommandLineProcessor;
00077   using Teuchos::RCP;
00078  
00079   bool success = true;
00080   bool result;
00081   
00082   //
00083   // (A) Setup and get basic MPI info
00084   //
00085   
00086   Teuchos::GlobalMPISession mpiSession(&argc,&argv);
00087 #ifdef HAVE_MPI
00088   MPI_Comm mpiComm = MPI_COMM_WORLD;
00089 #endif
00090   
00091   //
00092   // (B) Setup the output stream (do output only on root process!)
00093   //
00094 
00095   Teuchos::RCP<Teuchos::FancyOStream>
00096     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00097   
00098   try {
00099 
00100     //
00101     // (C) Read in commandline options
00102     //
00103 
00104 
00105     CommandLineProcessor  clp;
00106     clp.throwExceptions(false);
00107     clp.addOutputSetupOptions(true);
00108 
00109     int globalDim = 500;
00110     clp.setOption( "global-dim", &globalDim, "Global dimension of the linear system." );
00111 
00112     bool dumpAll = false;
00113     clp.setOption( "dump-all", "no-dump", &dumpAll, "Determines if quantities are dumped or not." );
00114 
00115     CommandLineProcessor::EParseCommandLineReturn
00116       parse_return = clp.parse(argc,argv);
00117     if (parse_return != CommandLineProcessor::PARSE_SUCCESSFUL)
00118       return parse_return;
00119 
00120     TEST_FOR_EXCEPTION( globalDim < 2, std::logic_error, "Error, globalDim=" << globalDim << " < 2 is not allowed!" );
00121 
00122     *out << "\n***\n*** Running power method example using Epetra implementation\n***\n" << std::scientific;
00123 
00124     //
00125     // (D) Setup the operator and run the power method!
00126     //
00127 
00128     //
00129     // (1) Setup the initial tridiagonal operator
00130     //
00131     //       [  2  -1             ]
00132     //       [ -1   2  -1         ]
00133     //  A =  [      .   .   .     ]
00134     //       [          -1  2  -1 ]
00135     //       [             -1   2 ]
00136     //
00137     *out << "\n(1) Constructing tridagonal Epetra matrix A of global dimension = " << globalDim << " ...\n";
00138     RCP<Epetra_Operator>
00139       A_epetra = createTridiagEpetraLinearOp(
00140         globalDim,
00141 #ifdef HAVE_MPI
00142         mpiComm,
00143 #endif
00144         1.0, true, *out
00145         );
00146     // Wrap in an Thyra::EpetraLinearOp object
00147     RCP<Thyra::LinearOpBase<double> >
00148       A = Thyra::nonconstEpetraLinearOp(A_epetra);
00149     //
00150     if (dumpAll) *out << "\nA =\n" << *A; // This works even in parallel!
00151     
00152     //
00153     // (2) Run the power method ANA
00154     //
00155     *out << "\n(2) Running the power method on matrix A ...\n";
00156     double  lambda      = 0.0;
00157     double  tolerance   = 1e-3;
00158     int     maxNumIters = 10*globalDim;
00159     result = sillyPowerMethod<double>(*A, maxNumIters, tolerance, outArg(lambda), *out);
00160     if(!result) success = false;
00161     *out << "\n  Estimate of dominate eigenvalue lambda = " << lambda << std::endl;
00162     
00163     //
00164     // (3) Increase dominance of first eigenvalue
00165     //
00166     *out << "\n(3) Scale the diagonal of A by a factor of 10 ...\n";
00167     scaleFirstDiagElement( 10.0, &*A );
00168     
00169     //
00170     // (4) Run the power method ANA again
00171     //
00172     *out << "\n(4) Running the power method again on matrix A ...\n";
00173     result = sillyPowerMethod<double>(*A, maxNumIters, tolerance, outArg(lambda), *out);
00174     if(!result) success = false;
00175     *out << "\n  Estimate of dominate eigenvalue lambda = " << lambda << std::endl;
00176     
00177   }
00178   TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success)
00179   
00180   if (success)
00181     *out << "\nCongratulations! All of the tests checked out!\n";
00182   else
00183     *out << "\nOh no! At least one of the tests failed!\n";
00184 
00185   return success ? 0 : 1;
00186 
00187 } // end main()
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines