sillyCgSolve_mpi.cpp

Click here for a more detailed discussion of this example program.

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 // 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 "ExampleTridiagSpmdLinearOp.hpp"
00030 #include "sillyCgSolve.hpp"
00031 #include "Thyra_VectorStdOps.hpp"
00032 #include "Thyra_TestingTools.hpp"
00033 #include "Thyra_LinearOpTester.hpp"
00034 #include "Teuchos_GlobalMPISession.hpp"
00035 #include "Teuchos_CommandLineProcessor.hpp"
00036 #include "Teuchos_DefaultComm.hpp"
00037 #include "Teuchos_VerboseObject.hpp"
00038 #include "Teuchos_Time.hpp"
00039 #include "Teuchos_StandardCatchMacros.hpp"
00040 
00041 //
00042 // This example program is meant to show how easy it is to create MPI Thyra
00043 // objects and use them with an ANA (CG in this case).
00044 //
00045 // This example uses a silly concrete tridiagonal matrix class called
00046 // ExampleTridiagSpmdLinearOp that demonstrates how to write such subclasses.
00047 //
00048 template<class Scalar>
00049 bool runCgSolveExample(
00050   const Teuchos::RCP<const Teuchos::Comm<Thyra::Index> > &comm
00051   ,const int                                                     procRank
00052   ,const int                                                     numProc
00053   ,const int                                                     localDim
00054   ,const Scalar                                                  diagScale
00055   ,const bool                                                    showAllTests
00056   ,const bool                                                    verbose
00057   ,const bool                                                    dumpAll
00058   ,const typename Teuchos::ScalarTraits<Scalar>::magnitudeType   tolerance
00059   ,const int                                                     maxNumIters
00060   )
00061 {
00062   using Teuchos::RCP;
00063   using Teuchos::rcp;
00064   using Teuchos::OSTab;
00065   typedef Teuchos::ScalarTraits<Scalar> ST;
00066   typedef typename ST::magnitudeType    ScalarMag;
00067   bool success = true;
00068   bool result;
00069   // ToDo: Get VerboseObjectBase to automatically setup for parallel
00070   Teuchos::RCP<Teuchos::FancyOStream>
00071     out = (verbose ? Teuchos::VerboseObjectBase::getDefaultOStream() : Teuchos::null);
00072   if(verbose)
00073     *out << "\n***\n*** Running silly CG solver using scalar type = \'" << ST::name() << "\' ...\n***\n";
00074   Teuchos::Time timer("");
00075   timer.start(true);
00076   //
00077   // (A) Setup a simple linear system with tridiagonal operator:
00078   //
00079   //       [  a*2   -1                ]
00080   //       [ -1    a*2  -1            ]
00081   //  A =  [         .   .    .       ]
00082   //       [            -1  a*2    -1 ]
00083   //       [                 -1   a*2 ]
00084   //
00085   // (A.1) Create the tridiagonal matrix operator
00086   if(verbose)
00087     *out << "\nConstructing tridiagonal matrix A of local dimension = " << localDim
00088          << " and diagonal multiplier = " << diagScale << " ...\n";
00089   const Thyra::Index
00090     lowerDim = ( procRank == 0         ? localDim - 1 : localDim ),
00091     upperDim = ( procRank == numProc-1 ? localDim - 1 : localDim );
00092   std::vector<Scalar> lower(lowerDim), diag(localDim), upper(upperDim);
00093   const Scalar one = ST::one(), diagTerm = Scalar(2)*diagScale*ST::one();
00094   int k = 0, kl = 0;
00095   if(procRank > 0) { lower[kl] = -one; ++kl; };  diag[k] = diagTerm; upper[k] = -one; // First local row
00096   for( k = 1; k < localDim - 1; ++k, ++kl ) {
00097     lower[kl] = -one; diag[k] = diagTerm; upper[k] = -one;                            // Middle local rows
00098   }
00099   lower[kl] = -one; diag[k] = diagTerm; if(procRank < numProc-1) upper[k] = -one;     // Last local row
00100   RCP<const Thyra::LinearOpBase<Scalar> >
00101     A = rcp(new ExampleTridiagSpmdLinearOp<Scalar>(comm,localDim,&lower[0],&diag[0],&upper[0]));
00102   if(verbose) *out << "\nGlobal dimension of A = " << A->domain()->dim() << std::endl;
00103   // (A.2) Testing the linear operator constructed linear operator
00104   if(verbose) *out << "\nTesting the constructed linear operator A ...\n";
00105   Thyra::LinearOpTester<Scalar> linearOpTester;
00106   linearOpTester.dump_all(dumpAll);
00107   linearOpTester.set_all_error_tol(tolerance);
00108   linearOpTester.set_all_warning_tol(ScalarMag(ScalarMag(1e-2)*tolerance));
00109   linearOpTester.show_all_tests(true);
00110   linearOpTester.check_adjoint(false);
00111   linearOpTester.check_for_symmetry(true);
00112   linearOpTester.show_all_tests(showAllTests);
00113   result = linearOpTester.check(*A,out.get());
00114   if(!result) success = false;
00115   // (A.3) Create RHS vector b and set to a random value
00116   RCP<Thyra::VectorBase<Scalar> > b = createMember(A->range());
00117   Thyra::seed_randomize<Scalar>(0);
00118   Thyra::randomize( Scalar(-ST::one()), Scalar(+ST::one()), &*b );
00119   // (A.4) Create LHS vector x and set to zero
00120   RCP<Thyra::VectorBase<Scalar> > x = createMember(A->domain());
00121   Thyra::assign( &*x, ST::zero() );
00122   //
00123   // (B) Solve the linear system with the silly CG solver
00124   //
00125   if(verbose) *out << "\nSolving the linear system with sillyCgSolve(...) ...\n";
00126   result = sillyCgSolve(*A,*b,maxNumIters,tolerance,&*x,OSTab(out).get());
00127   if(!result) success = false;
00128   //
00129   // (C) Check that the linear system was solved to the specified tolerance
00130   //
00131   RCP<Thyra::VectorBase<Scalar> > r = createMember(A->range());                     
00132   Thyra::assign(&*r,*b);                                               // r = b
00133   Thyra::apply(*A,Thyra::NOTRANS,*x,&*r,Scalar(-ST::one()),ST::one()); // r = -A*x + r
00134   const ScalarMag r_nrm = Thyra::norm(*r), b_nrm = Thyra::norm(*b);
00135   const ScalarMag rel_err = r_nrm/b_nrm, relaxTol = ScalarMag(10.0)*tolerance;
00136   result = rel_err <= relaxTol;
00137   if(!result) success = false;
00138   if(verbose) {
00139     *out << "\nChecking the residual ourselves ...\n";
00140     {
00141       OSTab tab(out);
00142       *out
00143         << "\n||b-A*x||/||b|| = "<<r_nrm<<"/"<<b_nrm<<" = "<<rel_err<<(result?" <= ":" > ")
00144         <<"10.0*tolerance = "<<relaxTol<<": "<<(result?"passed":"failed")<<std::endl;
00145     }
00146   }
00147   timer.stop();
00148   if(verbose) *out << "\nTotal time = " << timer.totalElapsedTime() << " sec\n";
00149 
00150   return success;
00151 } // end runCgSolveExample()
00152 
00153 //
00154 // Actual main driver program
00155 //
00156 int main(int argc, char *argv[])
00157 {
00158 
00159   using Teuchos::CommandLineProcessor;
00160 
00161   bool success = true;
00162   bool verbose = true;
00163   bool result;
00164 
00165   Teuchos::GlobalMPISession mpiSession(&argc,&argv);
00166   const int procRank = Teuchos::GlobalMPISession::getRank();
00167   const int numProc = Teuchos::GlobalMPISession::getNProc();
00168 
00169   const Teuchos::RCP<const Teuchos::Comm<Thyra::Index> >
00170     comm = Teuchos::DefaultComm<Thyra::Index>::getComm();
00171 
00172   Teuchos::RCP<Teuchos::FancyOStream>
00173     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00174 
00175   try {
00176 
00177     //
00178     // Read in command-line options
00179     //
00180 
00181     int    localDim    = 500;
00182     double diagScale   = 1.001;
00183     double tolerance   = 1e-4;
00184     bool   showAllTests = false;
00185     int    maxNumIters = 300;
00186     bool   dumpAll     = false;
00187 
00188     CommandLineProcessor  clp;
00189     clp.throwExceptions(false);
00190     clp.addOutputSetupOptions(true);
00191 
00192     clp.setOption( "verbose", "quiet", &verbose, "Determines if any output is printed or not." );
00193     clp.setOption( "local-dim", &localDim, "Local dimension of the linear system." );
00194     clp.setOption( "diag-scale", &diagScale, "Scaling of the diagonal to improve conditioning." );
00195     clp.setOption( "show-all-tests", "show-summary-only", &showAllTests, "Show all LinearOpTester tests or not" );
00196     clp.setOption( "tol", &tolerance, "Relative tolerance for linear system solve." );
00197     clp.setOption( "max-num-iters", &maxNumIters, "Maximum of CG iterations." );
00198     clp.setOption( "dump-all", "no-dump-all", &dumpAll, "Determines if vectors are printed or not." );
00199 
00200     CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
00201     if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
00202 
00203     TEST_FOR_EXCEPTION( localDim < 2, std::logic_error, "Error, localDim=" << localDim << " < 2 is not allowed!" );
00204 
00205     // Run using float
00206     result = runCgSolveExample<float>(comm,procRank,numProc,localDim,diagScale,showAllTests,verbose,dumpAll,tolerance,maxNumIters);
00207     if(!result) success = false;
00208 
00209     // Run using double
00210     result = runCgSolveExample<double>(comm,procRank,numProc,localDim,diagScale,showAllTests,verbose,dumpAll,tolerance,maxNumIters);
00211     if(!result) success = false;
00212 
00213 #if defined(HAVE_COMPLEX) && defined(HAVE_TEUCHOS_COMPLEX)
00214 
00215     // Run using std::complex<float>
00216     result = runCgSolveExample<std::complex<float> >(comm,procRank,numProc,localDim,diagScale,showAllTests,verbose,dumpAll,tolerance,maxNumIters);
00217     if(!result) success = false;
00218 
00219     // Run using std::complex<double>
00220     result = runCgSolveExample<std::complex<double> >(comm,procRank,numProc,localDim,diagScale,showAllTests,verbose,dumpAll,tolerance,maxNumIters);
00221     if(!result) success = false;
00222 
00223 #endif    
00224 
00225   }
00226   TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success)
00227 
00228   if( verbose && procRank==0 ) {
00229     if(success) *out << "\nAll of the tests seem to have run successfully!\n";
00230     else        *out << "\nOh no! at least one of the tests failed!\n"; 
00231   }
00232   
00233   return success ? 0 : 1;
00234 
00235 } // end main()

Generated on Wed May 12 21:42:27 2010 for Thyra Operator/Vector Support by  doxygen 1.4.7