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