Thyra Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
sillyCgSolve_epetra.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Thyra: Interfaces and Support for Abstract Numerical Algorithms
5// Copyright (2004) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov)
38//
39// ***********************************************************************
40// @HEADER
41
43#include "sillyCgSolve.hpp"
46#include "Thyra_DefaultSpmdVectorSpace.hpp"
51
52//
53// Main driver program for epetra implementation of CG.
54//
55int main(int argc, char *argv[])
56{
58 using Teuchos::RCP;
59
60 bool success = true;
61 bool verbose = true;
62 bool result;
63
64 //
65 // (A) Setup and get basic MPI info
66 //
67
68 Teuchos::GlobalMPISession mpiSession(&argc,&argv);
69#ifdef HAVE_MPI
70 MPI_Comm mpiComm = MPI_COMM_WORLD;
71#endif
72
73 //
74 // (B) Setup the output stream (do output only on root process!)
75 //
76
79
80 try {
81
82 //
83 // (C) Read in commandline options
84 //
85
86 int globalDim = 500;
87 double diagScale = 1.001;
88 bool useWithNonEpetraVectors = false;
89 double tolerance = 1e-4;
90 int maxNumIters = 300;
91
92 CommandLineProcessor clp(false); // Don't throw exceptions
93 clp.setOption( "verbose", "quiet", &verbose, "Determines if any output is printed or not." );
94 clp.setOption( "global-dim", &globalDim, "Global dimension of the linear system." );
95 clp.setOption( "diag-scale", &diagScale, "Scaling of the diagonal to improve conditioning." );
96 clp.setOption( "use-with-non-epetra-vectors", "use-with-epetra-vectors", &useWithNonEpetraVectors
97 , "Use non-epetra vectors with Epetra operator or not." );
98 clp.setOption( "tol", &tolerance, "Relative tolerance for linear system solve." );
99 clp.setOption( "max-num-iters", &maxNumIters, "Maximum of CG iterations." );
100
101 CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
102 if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
103
104 TEUCHOS_TEST_FOR_EXCEPTION( globalDim < 2, std::logic_error, "Error, globalDim=" << globalDim << " < 2 is not allowed!" );
105
106 if(verbose) *out << "\n***\n*** Running CG example using Epetra implementation\n***\n" << std::scientific;
107
108 //
109 // (D) Setup a simple linear system with tridagonal operator:
110 //
111 // [ a*2 -1 ]
112 // [ -1 a*2 -1 ]
113 // A = [ . . . ]
114 // [ -1 a*2 -1 ]
115 // [ -1 a*2 ]
116 //
117 // (D.1) Create the tridagonal Epetra matrix operator
118 if(verbose) *out << "\n(1) Constructing tridagonal Epetra matrix A_epetra of global dimension = " << globalDim << " ...\n";
119 RCP<Epetra_Operator>
121 globalDim
122#ifdef HAVE_MPI
123 ,mpiComm
124#endif
125 ,diagScale,verbose,*out
126 );
127 // (D.2) Wrap the Epetra_Opertor in a Thyra::EpetraLinearOp object
128 RCP<const Thyra::LinearOpBase<double> >
129 A = Thyra::epetraLinearOp(A_epetra);
130 // (D.3) Create RHS vector b and set to a random value
131 RCP<const Thyra::VectorSpaceBase<double> >
132 b_space = A->range();
133 // (D.4) Create the RHS vector b and initialize it to a random vector
134 RCP<Thyra::VectorBase<double> > b = createMember(b_space);
136 Thyra::randomize( -1.0, +1.0, b.ptr() );
137 // (D.5) Create LHS vector x and set to zero
138 RCP<Thyra::VectorBase<double> > x = createMember(A->domain());
139 Thyra::assign( x.ptr(), 0.0 );
140 //
141 // (E) Solve the linear system with the silly CG solver
142 //
143 result = sillyCgSolve(*A, *b, maxNumIters, tolerance, x.ptr(), *out);
144 if(!result) success = false;
145 //
146 // (F) Check that the linear system was solved to the specified tolerance
147 //
148 RCP<Thyra::VectorBase<double> > r = createMember(A->range());
149 Thyra::assign(r.ptr(),*b); // r = b
150 Thyra::apply(*A,Thyra::NOTRANS,*x,r.ptr(),-1.0,+1.0); // r = -A*x + r
151 const double r_nrm = Thyra::norm(*r), b_nrm =Thyra:: norm(*b);
152 const double rel_err = r_nrm/b_nrm, relaxTol = 2.0*tolerance;
153 result = rel_err <= relaxTol;
154 if(!result) success = false;
155 if(verbose)
156 *out
157 << "\n||b-A*x||/||b|| = "<<r_nrm<<"/"<<b_nrm<<" = "<<rel_err<<(result?" <= ":" > ")
158 <<"2.0*tolerance = "<<relaxTol<<": "<<(result?"passed":"failed")<<std::endl;
159
160 }
161 TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success)
162
163 if (verbose) {
164 if(success) *out << "\nCongratulations! All of the tests checked out!\n";
165 else *out << "\nOh no! At least one of the tests failed!\n";
166 }
167
168 return success ? 0 : 1;
169
170} // end main()
int main(int argc, char *argv[])
#define TEUCHOS_STANDARD_CATCH_STATEMENTS(VERBOSE, ERR_STREAM, SUCCESS_FLAG)
static RCP< FancyOStream > getDefaultOStream()
Teuchos::RCP< Epetra_Operator > createTridiagEpetraLinearOp(const int globalDim, const double diagScale, const bool verbose, std::ostream &out)
This function generates a tridiagonal linear operator using Epetra.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)