Anasazi Version of the Day
Loading...
Searching...
No Matches
AnasaziStatusTestSpecTrans.hpp
1// @HEADER
2// ***********************************************************************
3//
4// Anasazi: Block Eigensolvers Package
5// Copyright 2004 Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
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 Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40// @HEADER
41
42#ifndef ANASAZI_STATUS_TEST_SPECTRANS_HPP
43#define ANASAZI_STATUS_TEST_SPECTRANS_HPP
44
45#include "AnasaziTypes.hpp"
48
49using Teuchos::RCP;
50
51namespace Anasazi {
52namespace Experimental {
53
54 template<class ScalarType, class MV, class OP>
55 class StatusTestSpecTrans : public StatusTest<ScalarType,MV,OP> {
56
57 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
58 typedef MultiVecTraits<ScalarType,MV> MVT;
59 typedef OperatorTraits<ScalarType,MV,OP> OPT;
60
61 public:
62
63 // Constructor
64 StatusTestSpecTrans(MagnitudeType tol, int quorum = -1, ResType whichNorm = RES_2NORM, bool scaled = true, bool throwExceptionOnNan = true, const RCP<const OP> Mop = Teuchos::null);
65
66 // Destructor
67 virtual ~StatusTestSpecTrans() {};
68
69 // Check whether the test passed or failed
70 TestStatus checkStatus(Eigensolver<ScalarType,MV,OP> *solver);
71
72 // Return the result of the most recent checkStatus call
73 TestStatus getStatus() const { return state_; }
74
75 // Get the indices for the vectors that passed the test
76 std::vector<int> whichVecs() const { return ind_; }
77
78 // Get the number of vectors that passed the test
79 int howMany() const { return ind_.size(); }
80
81 void setQuorum (int quorum) {
82 state_ = Undefined;
83 quorum_ = quorum;
84 }
85
86 int getQuorum() const { return quorum_; }
87
88 void setTolerance(MagnitudeType tol)
89 {
90 state_ = Undefined;
91 tol_ = tol;
92 }
93
94 MagnitudeType getTolerance() const { return tol_; }
95
96 void setWhichNorm(ResType whichNorm)
97 {
98 state_ = Undefined;
99 whichNorm_ = whichNorm;
100 }
101
102 ResType getWhichNorm() const { return whichNorm_; }
103
104 void setScale(bool relscale)
105 {
106 state_ = Undefined;
107 scaled_ = relscale;
108 }
109
110 bool getScale() const { return scaled_; }
111
112 // Informs the status test that it should reset its internal configuration to the uninitialized state
113 void reset()
114 {
115 ind_.resize(0);
116 state_ = Undefined;
117 }
118
119 // Clears the results of the last status test
120 void clearStatus() { reset(); };
121
122 // Output formatted description of stopping test to output stream
123 std::ostream & print(std::ostream &os, int indent=0) const;
124
125 private:
126 TestStatus state_;
127 MagnitudeType tol_;
128 std::vector<int> ind_;
129 int quorum_;
130 bool scaled_;
131 enum ResType whichNorm_;
132 bool throwExceptionOnNaN_;
133 RCP<const OP> M_;
134
135 const MagnitudeType ONE;
136 };
137
138
139
140 template <class ScalarType, class MV, class OP>
141 StatusTestSpecTrans<ScalarType,MV,OP>::StatusTestSpecTrans(MagnitudeType tol, int quorum, ResType whichNorm, bool scaled, bool throwExceptionOnNaN, const RCP<const OP> Mop)
142 : state_(Undefined),
143 tol_(tol),
144 quorum_(quorum),
145 scaled_(scaled),
146 whichNorm_(whichNorm),
147 throwExceptionOnNaN_(throwExceptionOnNaN),
148 M_(Mop),
149 ONE(Teuchos::ScalarTraits<MagnitudeType>::one())
150 {}
151
152
153
154 template <class ScalarType, class MV, class OP>
155 TestStatus StatusTestSpecTrans<ScalarType,MV,OP>::checkStatus( Eigensolver<ScalarType,MV,OP>* solver )
156 {
157 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
158 typedef TraceMinBase<ScalarType,MV,OP> TS;
159
160 // Cast the eigensolver to a TraceMin solver
161 // Throw an exception if this fails
162 TS* tm_solver = dynamic_cast<TS*>(solver);
163 TEUCHOS_TEST_FOR_EXCEPTION(tm_solver == 0, std::invalid_argument, "The status test for spectral transformations currently only works for the trace minimization eigensolvers. Sorry!");
164
165 // Get the residual Ax-\lambda Bx, which is not computed when there's a spectral transformation
166 // TraceMin computes Bx-1/\lambda Ax
167 TraceMinBaseState<ScalarType,MV> state = tm_solver->getState();
168
169 size_t nvecs = state.ritzShifts->size();
170 std::vector<int> curind(nvecs);
171 for(size_t i=0; i<nvecs; i++)
172 curind[i] = i;
173
174 RCP<const MV> locKX, locMX, locX;
175 RCP<MV> R;
176 locX = MVT::CloneView(*state.X,curind);
177 if(state.KX != Teuchos::null)
178 locKX = MVT::CloneView(*state.KX,curind);
179 else
180 locKX = locX;
181 if(state.MX != Teuchos::null)
182 locMX = MVT::CloneView(*state.MX,curind);
183 else
184 locMX = locX;
185 R = MVT::CloneCopy(*locKX,curind);
186
187 std::vector<MagnitudeType> evals(nvecs);
188 for(size_t i=0; i<nvecs; i++)
189 evals[i] = ONE/(*state.T)[i];
190 MVT::MvScale(*R,evals);
191 MVT::MvAddMv(-ONE,*R,ONE,*locMX,*R);
192
193 // Compute the norms
194 std::vector<MagnitudeType> res(nvecs);
195 switch (whichNorm_) {
196 case RES_2NORM:
197 {
198 MVT::MvNorm(*R,res);
199 break;
200 }
201 case RES_ORTH:
202 {
203 RCP<MV> MR = MVT::Clone(*R,nvecs);
204 OPT::Apply(*M_,*R,*MR);
205 MVT::MvDot(*R,*MR,res);
206 for(size_t i=0; i<nvecs; i++)
207 res[i] = MT::squareroot(res[i]);
208 break;
209 }
210 case RITZRES_2NORM:
211 {
212 TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, "The trace minimization solvers do not define a Ritz residual. Please choose a different residual type");
213 break;
214 }
215 }
216
217 // if appropriate, scale the norms by the magnitude of the eigenvalue estimate
218 if(scaled_)
219 {
220 for(size_t i=0; i<nvecs; i++)
221 res[i] /= std::abs(evals[i]);
222 }
223
224 // test the norms
225 ind_.resize(0);
226 for(size_t i=0; i<nvecs; i++)
227 {
228 TEUCHOS_TEST_FOR_EXCEPTION( MT::isnaninf(res[i]), ResNormNaNError,
229 "StatusTestSpecTrans::checkStatus(): residual norm is nan or inf" );
230 if(res[i] < tol_)
231 ind_.push_back(i);
232 }
233 int have = ind_.size();
234 int need = (quorum_ == -1) ? nvecs : quorum_;
235 state_ = (have >= need) ? Passed : Failed;
236 return state_;
237 }
238
239
240
241 template <class ScalarType, class MV, class OP>
242 std::ostream& StatusTestSpecTrans<ScalarType,MV,OP>::print(std::ostream& os, int indent) const
243 {
244 std::string ind(indent,' ');
245 os << ind << "- StatusTestSpecTrans: ";
246 switch (state_) {
247 case Passed:
248 os << "Passed\n";
249 break;
250 case Failed:
251 os << "Failed\n";
252 break;
253 case Undefined:
254 os << "Undefined\n";
255 break;
256 }
257 os << ind << " (Tolerance, WhichNorm,Scaled,Quorum): "
258 << "(" << tol_;
259 switch (whichNorm_) {
260 case RES_ORTH:
261 os << ",RES_ORTH";
262 break;
263 case RES_2NORM:
264 os << ",RES_2NORM";
265 break;
266 case RITZRES_2NORM:
267 os << ",RITZRES_2NORM";
268 break;
269 }
270 os << "," << (scaled_ ? "true" : "false")
271 << "," << quorum_
272 << ")\n";
273
274 if (state_ != Undefined) {
275 os << ind << " Which vectors: ";
276 if (ind_.size() > 0) {
277 for(size_t i=0; i<ind_.size(); i++) os << ind_[i] << " ";
278 os << std::endl;
279 }
280 else
281 os << "[empty]\n";
282 }
283 return os;
284 }
285
286}} // end of namespace
287
288#endif
A status test for testing the norm of the eigenvectors residuals.
Abstract base class for trace minimization eigensolvers.
Types and exceptions used within Anasazi solvers and interfaces.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.
ResType
Enumerated type used to specify which residual norm used by residual norm status tests.
TestStatus
Enumerated type used to pass back information from a StatusTest.