sillyCgSolve_serial.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 "ExampleTridiagSerialLinearOp.hpp"
00030 #include "sillyCgSolve.hpp"
00031 #include "sillierCgSolve.hpp"
00032 #include "Thyra_DefaultScaledAdjointLinearOp.hpp"
00033 #include "Thyra_DefaultMultipliedLinearOp.hpp"
00034 #include "Thyra_VectorStdOps.hpp"
00035 #include "Thyra_TestingTools.hpp"
00036 #include "Thyra_LinearOpTester.hpp"
00037 #include "Teuchos_GlobalMPISession.hpp"
00038 #include "Teuchos_CommandLineProcessor.hpp"
00039 #include "Teuchos_VerboseObject.hpp"
00040 #include "Teuchos_Time.hpp"
00041 #include "Teuchos_StandardCatchMacros.hpp"
00042 
00043 
00044 //
00045 // This example program is meant to show how easy it is to create
00046 // serial Thyra objects and use them with an ANA (CG in this case).
00047 //
00048 // This example uses a silly concrete tridiagonal matrix class called
00049 // ExampleSpmdTridiagLinearOp that demonstrates how to write and use such
00050 // subclasses.
00051 //
00052 template<class Scalar>
00053 bool runCgSolveExample(
00054   const int dim,
00055   const Scalar diagScale,
00056   const bool symOp,
00057   const bool showAllTests,
00058   const typename Teuchos::ScalarTraits<Scalar>::magnitudeType tolerance,
00059   const int maxNumIters,
00060   const bool useSillierCg
00061   )
00062 {
00063 
00064   using Teuchos::as;
00065   using Teuchos::null;
00066   using Teuchos::RCP;
00067   using Teuchos::rcp;
00068   using Teuchos::OSTab;
00069   typedef Teuchos::ScalarTraits<Scalar> ST;
00070   using Thyra::multiply;
00071   using Thyra::scale;
00072   typedef typename ST::magnitudeType  ScalarMag;
00073   bool success = true;
00074   bool result;
00075   Teuchos::RCP<Teuchos::FancyOStream> out =
00076     Teuchos::VerboseObjectBase::getDefaultOStream();
00077   *out << "\n***\n*** Running silly CG solver using scalar type = \'"
00078        << ST::name() << "\' ...\n***\n";
00079   Teuchos::Time timer("");
00080   timer.start(true);
00081 
00082   //
00083   // (A) Setup a simple linear system with tridiagonal operator:
00084   //
00085   //       [   a*2   -1                         ]
00086   //       [ -r(1)  a*2       -1                ]
00087   //  A =  [          .        .        .       ]
00088   //       [             -r(n-2)      a*2    -1 ]
00089   //       [                      -r(n-1)   a*2 ]
00090   //
00091 
00092   // (A.1) Create the tridiagonal matrix operator
00093   *out << "\nConstructing tridiagonal matrix A of dimension = " << dim
00094        << " and diagonal multiplier = " << diagScale << " ...\n";
00095   std::vector<Scalar> lower(dim-1), diag(dim), upper(dim-1);
00096   const Scalar
00097     up = -ST::one(),
00098     diagTerm = as<Scalar>(2.0) * diagScale * ST::one(),
00099     low = -(symOp ? ST::one() : ST::random());
00100   int k = 0;
00101   // First row
00102   diag[k] = diagTerm; upper[k] = up;
00103   // Middle rows
00104   for( k = 1; k < dim - 1; ++k ) {
00105     lower[k-1] = low; diag[k] = diagTerm; upper[k] = up;
00106   }
00107   // Last row
00108   lower[k-1] = low; diag[k] = diagTerm;
00109   RCP<const Thyra::LinearOpBase<Scalar> > A =
00110     rcp(new ExampleTridiagSerialLinearOp<Scalar>(dim, &lower[0], &diag[0], &upper[0]));
00111 
00112   // (A.2) Testing the linear operator constructed linear operator
00113   *out << "\nTesting the constructed linear operator A ...\n";
00114   Thyra::LinearOpTester<Scalar> linearOpTester;
00115   linearOpTester.enable_all_tests(false);
00116   linearOpTester.check_linear_properties(true);
00117   linearOpTester.set_all_error_tol(tolerance);
00118   linearOpTester.set_all_warning_tol(1e-2*tolerance);
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   // (A.5) Create the final linear system
00133   if(!symOp) {
00134     *out << "\nSetting up normal equations for unsymmetric system A^H*(A*x-b) => new A*x = b ...\n";
00135     // A^H*A
00136     RCP<const Thyra::LinearOpBase<Scalar> > AtA = multiply(adjoint(A), A);
00137     // A^H*b
00138     RCP<Thyra::VectorBase<Scalar> > nb = createMember(AtA->range());
00139     Thyra::apply<Scalar>(*A, Thyra::CONJTRANS, *b, nb.ptr());
00140     A = AtA;
00141     b = nb;
00142   }
00143 
00144   // (A.6) Testing the linear operator used with the solve
00145   *out << "\nTesting the linear operator used with the solve ...\n";
00146   linearOpTester.check_for_symmetry(true);
00147   result = linearOpTester.check(*A, out.ptr());
00148   if(!result) success = false;
00149 
00150   //
00151   // (B) Solve the linear system with the silly CG solver
00152   //
00153   *out << "\nSolving the linear system with sillyCgSolve(...) ...\n";
00154   if(useSillierCg) {
00155     OSTab tab(out);
00156     result = sillierCgSolve(*A, *b, maxNumIters, tolerance, x.ptr(), *out);
00157   }
00158   else {
00159     OSTab tab(out);
00160     result = sillyCgSolve(*A, *b, maxNumIters, tolerance, x.ptr(), *out);
00161   }
00162   if(!result) success = false;
00163 
00164   //
00165   // (C) Check that the linear system was solved to the specified tolerance
00166   //
00167   RCP<Thyra::VectorBase<Scalar> > r = createMember(A->range());                     
00168   // r = b
00169   Thyra::V_V(r.ptr(), *b);
00170    // r = -A*x + r
00171   Thyra::apply<Scalar>(*A, Thyra::NOTRANS, *x, r.ptr(), -ST::one(), ST::one());
00172   const ScalarMag r_nrm = Thyra::norm(*r), b_nrm = Thyra::norm(*b);
00173   const ScalarMag rel_err = r_nrm/b_nrm, relaxTol = 10.0*tolerance;
00174   result = rel_err <= relaxTol;
00175   if(!result) success = false;
00176   *out << "\nChecking the residual ourselves ...\n";
00177   {
00178     OSTab tab(out);
00179     *out
00180       << "\n||b-A*x||/||b|| = "<<r_nrm<<"/"<<b_nrm<<" = "<<rel_err<<(result?" <= ":" > ")
00181       <<"10.0*tolerance = "<<relaxTol<<": "<<(result?"passed":"failed")<<std::endl;
00182   }
00183   timer.stop();
00184   *out << "\nTotal time = " << timer.totalElapsedTime() << " sec\n";
00185   
00186   return success;
00187 
00188 } // end runCgSolveExample()
00189 
00190 
00191 //
00192 // Actual main driver program
00193 //
00194 int main(int argc, char *argv[])
00195 {
00196   
00197   using Teuchos::CommandLineProcessor;
00198 
00199   bool success = true;
00200   bool result;
00201 
00202   Teuchos::GlobalMPISession mpiSession(&argc,&argv);
00203   // Above, I only add this to allow this test to be run with an MPI build but
00204   // this test does not use MPI in any way.  This is needed with the MPICH
00205   // 1.2.7 used on by SIERRA.
00206 
00207   Teuchos::RCP<Teuchos::FancyOStream>
00208     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00209 
00210   try {
00211 
00212     //
00213     // Read in command-line options
00214     //
00215 
00216 
00217     CommandLineProcessor  clp;
00218     clp.throwExceptions(false);
00219     clp.addOutputSetupOptions(true);
00220 
00221     int dim = 500;
00222     clp.setOption( "dim", &dim,
00223       "Dimension of the linear system." );
00224 
00225     double diagScale = 1.001;
00226     clp.setOption( "diag-scale", &diagScale,
00227       "Scaling of the diagonal to improve conditioning." );
00228 
00229     bool symOp = true;
00230     clp.setOption( "sym-op", "unsym-op", &symOp,
00231       "Determines if the operator is symmetric or not." );
00232 
00233     bool showAllTests = false;
00234     clp.setOption( "show-all-tests", "show-summary-only", &showAllTests,
00235       "Show all LinearOpTester tests or not" );
00236 
00237     double tolerance = 1e-4;
00238     clp.setOption( "tol", &tolerance,
00239       "Relative tolerance for linear system solve." );
00240 
00241     int maxNumIters = 300;
00242     clp.setOption( "max-num-iters", &maxNumIters,
00243       "Maximum of CG iterations." );
00244 
00245     bool useSillierCg = false;
00246     clp.setOption( "use-sillier-cg", "use-silly-cg", &useSillierCg,
00247       "Use the handle-based sillerCgSolve() function or the nonhandle-based"
00248       " sillyCgSolve() function");
00249 
00250     CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
00251     if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
00252 
00253     TEST_FOR_EXCEPTION( dim < 2, std::logic_error, "Error, dim=" << dim << " < 2 is not allowed!" );
00254 
00255 #if defined(HAVE_THYRA_FLOAT)
00256     result = runCgSolveExample<float>(dim, diagScale, symOp, showAllTests,
00257       tolerance, maxNumIters, useSillierCg);
00258     if(!result) success = false;
00259 #endif
00260 
00261     result = runCgSolveExample<double>(dim, diagScale, symOp, showAllTests,
00262       tolerance, maxNumIters, useSillierCg);
00263     if(!result) success = false;
00264 
00265 #ifdef HAVE_THYRA_COMPLEX
00266 
00267 #if defined(HAVE_THYRA_FLOAT)
00268     result = runCgSolveExample<std::complex<float> >(dim, diagScale, symOp, showAllTests,
00269       tolerance, maxNumIters, useSillierCg);
00270     if(!result) success = false;
00271 #endif
00272 
00273     result = runCgSolveExample<std::complex<double> >(dim, diagScale, symOp, showAllTests,
00274       tolerance, maxNumIters, useSillierCg);
00275     if(!result) success = false;
00276 
00277 #endif // HAVE_THYRA_COMPLEX
00278 
00279 #ifdef HAVE_TEUCHOS_GNU_MP
00280 
00281     result = runCgSolveExample<mpf_class>(dim, diagScale, symOp, showAllTests, 
00282       tolerance, maxNumIters, useSillierCg);
00283     if(!result) success = false;
00284 
00285 #ifdef HAVE_THYRA_COMPLEX
00286 
00287     //result = runCgSolveExample<std::complex<mpf_class> >(dim, mpf_class(diagScale), symOp,
00288     showAllTests, mpf_class(tolerance), maxNumIters);
00289     //if(!result) success = false;
00290     //The above commented-out code throws a floating-point exception?
00291 
00292 #endif // HAVE_THYRA_COMPLEX
00293 
00294 #endif // HAVE_TEUCHOS_GNU_MP
00295 
00296   }
00297   TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success)
00298 
00299     if(success)
00300       *out << "\nCongratulations! All of the tests checked out!\n";
00301     else
00302       *out << "\nOh no! At least one of the tests failed!\n";
00303   
00304   return success ? 0 : 1;
00305 
00306 } // end main()

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