Thyra Package Browser (Single Doxygen Collection) Version of the Day
sillyCgSolve_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 "createTridiagEpetraLinearOp.hpp"
00043 #include "sillyCgSolve.hpp"
00044 #include "Thyra_EpetraThyraWrappers.hpp"
00045 #include "Thyra_EpetraLinearOp.hpp"
00046 #include "Thyra_DefaultSpmdVectorSpace.hpp"
00047 #include "Teuchos_GlobalMPISession.hpp"
00048 #include "Teuchos_CommandLineProcessor.hpp"
00049 #include "Teuchos_VerboseObject.hpp"
00050 #include "Teuchos_StandardCatchMacros.hpp"
00051 
00052 //
00053 // Main driver program for epetra implementation of CG.
00054 //
00055 int main(int argc, char *argv[])
00056 {
00057   using Teuchos::CommandLineProcessor;
00058   using Teuchos::RCP;
00059  
00060   bool success = true;
00061   bool verbose = true;
00062   bool result;
00063   
00064   //
00065   // (A) Setup and get basic MPI info
00066   //
00067 
00068   Teuchos::GlobalMPISession mpiSession(&argc,&argv);
00069 #ifdef HAVE_MPI
00070   MPI_Comm mpiComm = MPI_COMM_WORLD;
00071 #endif
00072 
00073   //
00074   // (B) Setup the output stream (do output only on root process!)
00075   //
00076   
00077   Teuchos::RCP<Teuchos::FancyOStream>
00078     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00079   
00080   try {
00081 
00082     //
00083     // (C) Read in commandline options
00084     //
00085 
00086     int    globalDim                  = 500;
00087     double diagScale                  = 1.001;
00088     bool   useWithNonEpetraVectors    = false;
00089     double tolerance                  = 1e-4;
00090     int    maxNumIters                = 300;
00091 
00092     CommandLineProcessor  clp(false); // Don't throw exceptions
00093     clp.setOption( "verbose", "quiet", &verbose, "Determines if any output is printed or not." );
00094     clp.setOption( "global-dim", &globalDim, "Global dimension of the linear system." );
00095     clp.setOption( "diag-scale", &diagScale, "Scaling of the diagonal to improve conditioning." );
00096     clp.setOption( "use-with-non-epetra-vectors", "use-with-epetra-vectors", &useWithNonEpetraVectors
00097                    , "Use non-epetra vectors with Epetra operator or not." );
00098     clp.setOption( "tol", &tolerance, "Relative tolerance for linear system solve." );
00099     clp.setOption( "max-num-iters", &maxNumIters, "Maximum of CG iterations." );
00100 
00101     CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
00102     if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
00103 
00104     TEUCHOS_TEST_FOR_EXCEPTION( globalDim < 2, std::logic_error, "Error, globalDim=" << globalDim << " < 2 is not allowed!" );
00105 
00106     if(verbose) *out << "\n***\n*** Running CG example using Epetra implementation\n***\n" << std::scientific;
00107 
00108     //
00109     // (D) Setup a simple linear system with tridagonal operator:
00110     //
00111     //       [  a*2   -1                ]
00112     //       [ -1    a*2  -1            ]
00113     //  A =  [         .   .    .       ]
00114     //       [            -1  a*2    -1 ]
00115     //       [                 -1   a*2 ]
00116     //
00117     // (D.1) Create the tridagonal Epetra matrix operator
00118     if(verbose) *out << "\n(1) Constructing tridagonal Epetra matrix A_epetra of global dimension = " << globalDim << " ...\n";
00119     RCP<Epetra_Operator>
00120       A_epetra = createTridiagEpetraLinearOp(
00121         globalDim
00122 #ifdef HAVE_MPI
00123         ,mpiComm
00124 #endif
00125         ,diagScale,verbose,*out
00126         );
00127     // (D.2) Wrap the Epetra_Opertor in a Thyra::EpetraLinearOp object
00128     RCP<const Thyra::LinearOpBase<double> >
00129       A = Thyra::epetraLinearOp(A_epetra);
00130     // (D.3) Create RHS vector b and set to a random value
00131     RCP<const Thyra::VectorSpaceBase<double> >
00132       b_space = A->range();
00133     // (D.4) Create the RHS vector b and initialize it to a random vector
00134     RCP<Thyra::VectorBase<double> > b = createMember(b_space);
00135     Thyra::seed_randomize<double>(0);
00136     Thyra::randomize( -1.0, +1.0, b.ptr() );
00137     // (D.5) Create LHS vector x and set to zero
00138     RCP<Thyra::VectorBase<double> > x = createMember(A->domain());
00139     Thyra::assign( x.ptr(), 0.0 );
00140     //
00141     // (E) Solve the linear system with the silly CG solver
00142     //
00143     result = sillyCgSolve(*A, *b, maxNumIters, tolerance, x.ptr(), *out);
00144     if(!result) success = false;
00145     //
00146     // (F) Check that the linear system was solved to the specified tolerance
00147     //
00148     RCP<Thyra::VectorBase<double> > r = createMember(A->range());                     
00149     Thyra::assign(r.ptr(),*b);                                        // r = b
00150     Thyra::apply(*A,Thyra::NOTRANS,*x,r.ptr(),-1.0,+1.0);             // r = -A*x + r
00151     const double r_nrm = Thyra::norm(*r), b_nrm =Thyra:: norm(*b);
00152     const double rel_err = r_nrm/b_nrm, relaxTol = 2.0*tolerance;
00153     result = rel_err <= relaxTol;
00154     if(!result) success = false;
00155     if(verbose)
00156       *out
00157         << "\n||b-A*x||/||b|| = "<<r_nrm<<"/"<<b_nrm<<" = "<<rel_err<<(result?" <= ":" > ")
00158         <<"2.0*tolerance = "<<relaxTol<<": "<<(result?"passed":"failed")<<std::endl;
00159     
00160   }
00161   TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success)
00162   
00163   if (verbose) {
00164     if(success)  *out << "\nCongratulations! All of the tests checked out!\n";
00165     else         *out << "\nOh no! At least one of the tests failed!\n";
00166   }
00167 
00168   return success ? 0 : 1;
00169 
00170 } // end main()
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines