Thyra Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
sillyPowerMethod_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
42#include "sillyPowerMethod.hpp"
44#include "Thyra_TestingTools.hpp"
51#include "Teuchos_dyn_cast.hpp"
52#include "Epetra_CrsMatrix.h"
53#include "Epetra_Map.h"
54
55//
56// Increase the first diagonal element of your tridiagonal matrix.
57//
58void scaleFirstDiagElement( const double diagScale, Thyra::LinearOpBase<double> *A )
59{
60 using Teuchos::RCP;
62 // (A) Get at the underlying Epetra_Operator object that the EpetraLinearOp
63 // object directly maintains.
64 const RCP<Epetra_Operator> epetra_op = Thyra::get_Epetra_Operator(*A);
65 // (B) Perform a dynamic cast to Epetra_CrsMatrix.
66 // Note, the dyn_cast<>() template function will throw std::bad_cast
67 // with a nice error message if the cast fails!
68 Epetra_CrsMatrix &crsMatrix = Teuchos::dyn_cast<Epetra_CrsMatrix>(*epetra_op);
69 // (C) Query if the first row of the matrix is on this processor and if
70 // it is get it and scale it.
71 if(crsMatrix.MyGlobalRow(0)) {
72 const int numRowNz = crsMatrix.NumGlobalEntries(0);
73 TEUCHOS_TEST_FOR_EXCEPT( numRowNz != 2 );
74 int returndNumRowNz; double rowValues[2]; int rowIndexes[2];
75 crsMatrix.ExtractGlobalRowCopy( 0, numRowNz, returndNumRowNz, &rowValues[0], &rowIndexes[0] );
76 TEUCHOS_TEST_FOR_EXCEPT( returndNumRowNz != 2 );
77 for( int k = 0; k < numRowNz; ++k) if (rowIndexes[k] == 0) rowValues[k] *= diagScale;
78 crsMatrix.ReplaceGlobalValues( 0, numRowNz, rowValues, rowIndexes );
79 }
80} // end scaleFirstDiagElement()
81
82//
83// Main driver program for epetra implementation of the power method.
84//
85int main(int argc, char *argv[])
86{
87
88 using Teuchos::outArg;
90 using Teuchos::RCP;
91
92 bool success = true;
93 bool result;
94
95 //
96 // (A) Setup and get basic MPI info
97 //
98
99 Teuchos::GlobalMPISession mpiSession(&argc,&argv);
100#ifdef HAVE_MPI
101 MPI_Comm mpiComm = MPI_COMM_WORLD;
102#endif
103
104 //
105 // (B) Setup the output stream (do output only on root process!)
106 //
107
110
111 try {
112
113 //
114 // (C) Read in commandline options
115 //
116
117
118 CommandLineProcessor clp;
119 clp.throwExceptions(false);
120 clp.addOutputSetupOptions(true);
121
122 int globalDim = 500;
123 clp.setOption( "global-dim", &globalDim, "Global dimension of the linear system." );
124
125 bool dumpAll = false;
126 clp.setOption( "dump-all", "no-dump", &dumpAll, "Determines if quantities are dumped or not." );
127
128 CommandLineProcessor::EParseCommandLineReturn
129 parse_return = clp.parse(argc,argv);
130 if (parse_return != CommandLineProcessor::PARSE_SUCCESSFUL)
131 return parse_return;
132
133 TEUCHOS_TEST_FOR_EXCEPTION( globalDim < 2, std::logic_error, "Error, globalDim=" << globalDim << " < 2 is not allowed!" );
134
135 *out << "\n***\n*** Running power method example using Epetra implementation\n***\n" << std::scientific;
136
137 //
138 // (D) Setup the operator and run the power method!
139 //
140
141 //
142 // (1) Setup the initial tridiagonal operator
143 //
144 // [ 2 -1 ]
145 // [ -1 2 -1 ]
146 // A = [ . . . ]
147 // [ -1 2 -1 ]
148 // [ -1 2 ]
149 //
150 *out << "\n(1) Constructing tridagonal Epetra matrix A of global dimension = " << globalDim << " ...\n";
151 RCP<Epetra_Operator>
153 globalDim,
154#ifdef HAVE_MPI
155 mpiComm,
156#endif
157 1.0, true, *out
158 );
159 // Wrap in an Thyra::EpetraLinearOp object
160 RCP<Thyra::LinearOpBase<double> >
161 A = Thyra::nonconstEpetraLinearOp(A_epetra);
162 //
163 if (dumpAll) *out << "\nA =\n" << *A; // This works even in parallel!
164
165 //
166 // (2) Run the power method ANA
167 //
168 *out << "\n(2) Running the power method on matrix A ...\n";
169 double lambda = 0.0;
170 double tolerance = 1e-3;
171 int maxNumIters = 10*globalDim;
172 result = sillyPowerMethod<double>(*A, maxNumIters, tolerance, outArg(lambda), *out);
173 if(!result) success = false;
174 *out << "\n Estimate of dominate eigenvalue lambda = " << lambda << std::endl;
175
176 //
177 // (3) Increase dominance of first eigenvalue
178 //
179 *out << "\n(3) Scale the diagonal of A by a factor of 10 ...\n";
180 scaleFirstDiagElement( 10.0, &*A );
181
182 //
183 // (4) Run the power method ANA again
184 //
185 *out << "\n(4) Running the power method again on matrix A ...\n";
186 result = sillyPowerMethod<double>(*A, maxNumIters, tolerance, outArg(lambda), *out);
187 if(!result) success = false;
188 *out << "\n Estimate of dominate eigenvalue lambda = " << lambda << std::endl;
189
190 }
191 TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success)
192
193 if (success)
194 *out << "\nCongratulations! All of the tests checked out!\n";
195 else
196 *out << "\nOh no! At least one of the tests failed!\n";
197
198 return success ? 0 : 1;
199
200} // 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_EXCEPT(throw_exception_test)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Teuchos::RCP< Epetra_Operator > get_Epetra_Operator(LinearOpBase< double > &op)
Full specialization for Scalar=double.
void scaleFirstDiagElement(const double diagScale, Thyra::LinearOpBase< double > *A)