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

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