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

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